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Abstract 


Numerical simulations of decaying homogeneous isotropic turbulence are performed with both low- 
order and high-order spatial discretization schemes. The turbulent Mach and Reynolds numbers 
for the simulations are 0.2 and 250, respectively. For the low-order schemes we use either second- 
order central or third-order upwind biased differencing. For higher order approximations we ap- 
ply weighted essentially non-oscillatory (WENO) schemes, both with linear and nonlinear weights. 
There are two objectives in this preliminary effort to investigate possible schemes for large eddy 
simulation (LES). One is to explore the capability of a widely used low-order computational fluid 
dynamics (CFD) code to perform LES computations. The other is to determine the effect of higher 
order accuracy (fifth, seventh, and ninth order) achieved with high-order upwind biased WENO- 
based schemes. Turbulence statistics, such as kinetic energy, dissipation, and skewness, along with 
the energy spectra from simulations of the decaying turbulence problem are used to assess and 
compare the various numerical schemes. In addition, results from the best performing schemes are 
compared with those from a spectral scheme. The effects of grid density, ranging from 32 cubed to 
192 cubed, on the computations are also examined. The fifth-order WENO-based scheme is found 
to be too dissipative, especially on the coarser grids. However, with the seventh-order and ninth- 
order WENO-based schemes we observe a significant improvement in accuracy relative to the lower 
order LES schemes, as revealed by the computed peak in the energy dissipation and by the energy 
spectrum. 


Nomenclature 

A = Jacobian matrix 

ADM = Approximate Deconvolution Model 
CFD = Computational Fluid Dynamics 
DIA = Direct Interaction Approximation 
DNS = Direct Numerical Simulation 
DSM = Differential Stress Model 
E = specific total energy 

E(k) = energy as a function of wave number in Fourier space 
ENO = Essentially Non-Oscillatory (numerical scheme) 

FCT = Flux-Corrected Transport (numerical scheme) 

G = LES convolution kernel 

Hi = sum of heat flux plus work done by stresses 

L = length scale 

LED = Local Extremum Diminishing 

LES = Large Eddy Simulation 

ILES = Implicit LES 

M = Mach number 

MILES = Monotone Implicit LES 

MUSCL = Monotone Upstream-centered Schemes for Conservation Laws (numerical scheme) 
PPM = Piecewise Parabolic Method (numerical scheme) 

Pr = Prandtl number 

Qi = subgrid scale heat flux 

R = gas constant 

Re = Reynolds number 

Re\ = Taylor microscale Reynolds number 

RK = Runge-Kutta 

S_={2S i jS ij ) 1 / 2 

Sij = strain rate tensor 
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Sk = velocity derivative skewness factor 

SGS = Sub-Grid Scale 

SSM = Scale-Similarity Model 

T = temperature 

Tj j = total stress tensor 

TVD = Total Variation Diminishing (numerical scheme) 

U = velocity 
Ui = state vector 

WENO = Weighted ENO (numerical scheme) 

WENO-L = Linear version of WENO (numerical scheme) 

c p = specific heat constant 

c s = constant in Smagorinsky SGS model 

c v = specific heat at constant volume 

e = internal energy per unit volume 

f = represents general variable; also represents numerical flux 
k = filtered stress kinetic energy; also used as an index counter 
p = pressure 
t = time 
t = time 

t* = time, nondimensionalized by U re f and L re f 
Ui = velocity component u, v, w 
x, = Cartesian component x, y, z 

Greek : 

A = LES filter width; also represents a spatial difference 
A = LES grid filter 

0 k = smoothness measure of WENO scheme 
Sij = Kronecker delta 
e = turbulent dissipation rate 

e w = parameter to prevent division by zero in WENO weights 
7 = ratio of specific heats 
// = Kolmogorov length scale 
k = wave number in Fourier space 
n' = coefficient of thermal conductivity 
k = variable in the MUSCL scheme 
A = Taylor microscale 
p = dynamic viscosity 
v = kinematic viscosity 
p = density 

Oij = laminar part of the stress tensor 
r = Kolmogorov time scale 
Tij = subgrid scale stress 

X = ratio of compressible kinetic energy to total kinetic energy 
u>k = weights of WENO scheme 
( t> = represents general space-time variable 

Subscripts, Superscripts, and Special Symbols : 

overbar = resolved part; also spatially-filtered variable 

tilde = Favre-averaged variable 

angled brackets = average over all cells in the domain 

C = compressible 

I = incompressible 
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L = left 

R = right 

re/ = reference 

rms = root-mean square 

t = turbulent 

x, y, z = in x, y , z directions 

0 = initial 

V = mathematical “del” differential operator 

1 Introduction 

1.1 General Description of DNS and LES 

Advancements in computer power, both processing speed and memory, in recent years and the 
need for improved representation of the effects of turbulence in practical engineering computations 
has encouraged increased interest in direct numerical simulation (DNS) and large eddy simulation 
(LES). In DNS the full Navier-Stokes equations are solved, and all the relevant scales, down to 
the Kolmogorov scale, of a physical problem are resolved. Unfortunately, even with the current 
computer capability, the computational requirements for DNS are too demanding for high Reynolds 
number (Re) problems since the number of grid points required is proportional to /ri: 9/1 . Further- 
more, as indicated by Piomelli [1], the computational cost scales like Re 3 . For these reasons DNS 
has generally been applied to flow problems with relatively low Re and simple geometries. Instead 
of resolving all the relevant scales of the turbulent motion, LES relies upon the assumption that the 
smaller scales of the turbulent motion have a more universal behavior that allows their effect (due 
to energy dissipation) on the larger scales to be represented by appropriate modeling. By modeling 
the physics of the smallest scales there is a significant reduction in the computational requirements 
relative to DNS, allowing the possibility to compute high Re turbulent flows over complex geo- 
metric configurations. There are two general approaches to LES. One involves the filtering of the 
Navier-Stokes equations in order to eliminate scales smaller than a prescribed value. The filtering 
introduces additional stresses and heat fluxes that are determined by the unresolved small scale (sub- 
grid scale) eddies. This approach relies on direct (explicit) modeling of the additional terms, and it 
is often called conventional or explicit LES. The second LES formulation depends on the numerical 
dissipation associated with the spatial discretization of the flow equations to emulate the physics of 
the smallest scales. This approach is called integrated or implicit LES (ILES). 

1.2 Conventional LES Methods 

In conventional LES the filtering to separate the large energy carrying eddies from the unresolved 
small scales is designed to remove the scales equal to or smaller than some filter width. This filter 
width is frequently represented by A, and it defines the characteristic resolution length. For the LES 
A //, with y = 0(Re~ 3 ' 4 ) denoting the Kolmogorov length scale, whereas for DNS there is no 
“filter” so it is equivalent to LES with Amy. The ultimate success in the LES computation depends 
upon how well the effects of the filtered scales on the larger resolved scales are represented. The ex- 
plicit subgrid-scale (SGS) models of conventional LES have two general classifications: functional 
models and structural models. There are also mixtures of these two types (e.g., functional models 
with insufficient dissipation of the small scales also require some type of additional dissipation). 

Functional models represent the subgrid terms (i.e., stress tensor and heat flux vector) indirectly 
by modeling the effect of the dissipative action of the small scales on the larger scales, whereas 
structural models represent the subgrid terms directly by approximating the unfiltered flow variables. 
One example of a functional model of the subgrid scale stress tensor is the Smagorinsky model [2], 
With this algebraic eddy viscosity model the SGS stress tensor is a function of the mean strain rate 


3 



tensor and the square of the filter width, which is usually taken as the mesh cell size. There is also a 
constant that must be specified. An improvement of this methodology is the dynamic Smagorinsky 
model [3], which replaces the square of the constant by a coefficient that is dynamically adjusted 
according to information available at the smallest resolved scales. In general, eddy-viscosity models 
can reproduce SGS dissipation but not the subgrid forces in the filtered momentum equations. 

In the case of the structural models the objective is to reconstruct some of the SGS information 
directly by generating a new solution field from the filtered one, resulting in an improved approxi- 
mation of the unfiltered field. Structural models are generally superior to functional models because 
they do not have the isotropic assumption that most functional models have. Some examples of 
structural models are the scale-similarity model (SSM) and mixed model of Bardina et al. [4], the 
differential stress equation model (DSM) of Deardotff [5], and the approximate deconvolution model 
(ADM) of Stolz and Adams [6,7]. In essence, the scale similarity hypothesis [4] assumes that the 
effect of the smallest resolved scales on the unresolved scales is similar to the effect of the largest 
resolved scales on the smallest resolved scales. Based on this similarity assumption, the SSM ap- 
proximates the SGS stress tensor using the resolved velocity field and a second application of the 
spatial filter. Due to the filtering process this model allows for backscatter (i.e., energy transfer 
from the subgrid scales to the resolved scales). The DSM [5] uses a transport equation to model 
the subgrid stress tensor. For the ADM [6, 7] the unfiltered velocity is estimated by an approximate 
deconvolution operation, which uses an approximation of the inverse filter kernel. Like the SSM, the 
ADM does not have adequate dissipation of energy. In order to represent the energy transfer from 
the resolved to unresolved scales, an additional relaxation is included in the discrete equations. This 
is equivalent to an explicit filtering of the resolved velocity field at each time step in the evolution 
of the solution. As pointed out by Grinstein and Fureby [8], a principal difficulty with explicit SGS 
models is the need to separate the effects of explicit filtering and SGS reconstruction models from 
those due to discretization. 

1.3 ILES Methods 

In general, one could classify any LES method that uses numerical dissipation rather than an explicit 
model to represent the effect of the subgrid scales as an ILES method. According to Sagaut [9] any 
ILES approach is implicitly based on the following hypothesis: The action of sub grid scales on the 
resolved scales is equivalent to a strictly dissipative action. The numerical dissipative terms are 
produced when approximating the convective terms in the fluid dynamic equations. These terms 
are usually either introduced implicitly through an upwind scheme framework or explicitly when 
using a central difference scheme. The numerical dissipation plays several roles. One role is to 
ensure numerical stability. In the linear sense, assuming well-posedness and consistency, stability 
is a necessary and sufficient condition for convergence (Lax equivalence theorem, see Richtmyer 
and Morton [10]). It is not a sufficient condition to guarantee physically realizable solutions. To 
ensure a unique and physically meaningful solution the dissipation must also allow satisfaction of 
an entropy condition. Another role for dissipation is to prevent oscillations in the solution when 
extrema occur due to physical phenomena such as shock waves and contact discontinuities. There is, 
in general, a reduction in the approximation order (e.g., second-order schemes become first order) for 
the convective terms at discontinuities in order to enable discontinuity capturing without oscillations. 
Frequently, the distinction between upwind and central difference methods depends primarily on 
the nature of the limiting (or switching) that is imposed to change the order of approximation at 
discontinuities. 

A subclass of ILES algorithms is the monotone integrated LES (MILES) schemes. To understand 
what methods can be placed in this category it is appropriate to provide some historical background 
and clarification. The MILES schemes are defined to be at least second-order accurate in smooth 
regions of the flow field. The classification of MILES was created in 1990 by Boris [1 1] in reference 
to the flux corrected transport (FCT) scheme introduced by Boris [12] in 1971 and generalized by 
Boris and Book [13] -[15], It stems from the fact that the FCT scheme has a predictor-corrector 
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structure with the first approximate solution determined with a first-order monotone scheme. This 
initial approximation is then corrected by removing the large dissipative error with a limited anti- 
diffusion. Through the limiting process monotonicity and positivity conditions are imposed. The 
monotonicity constraint, which is the underlying FCT principle given in [12] and [13], ensures 
that no new extrema are created and that there is no growth of existing extrema. The positivity 
condition prevents positive definite physical quantities (e.g., density) from becoming negative. The 
resulting algorithm is at least second-order accurate in smooth regions, and it satisfies monotonicity- 
preserving and positivity-preserving properties. In addition, the algorithm is designed to satisfy the 
physical properties of global conservation, causality, and locality. Causality simply means that a 
fluid passing between two locations in a discrete domain must pass through all the cells in the flow 
path between those locations. Since the FCT scheme is conservative, causality must be satisfied. 
Locality requires the approximations to derivatives to be compact. 

Boris classifies the FCT scheme as a monotone fluid dynamics algorithm. In [16] Boris lists 
other examples of monotone algorithms: monotone upstream-centered scheme for conservation 
laws (MUSCL) of van Leer [17] -[21]; piecewise parabolic method (PPM) of Woodward and 
Colella [22, 23]; total variation diminishing (TVD) schemes of Harten [24, 25], (TVD) schemes 
of Harten [24,25], Sweby [26], Yee et al. [27]; and second-order Godunov methods of Colella [28] 
and Colella and Glaz [29]. Monotone schemes are monotonicity preserving and guarantee a mono- 
tone solution if the scheme converges. In addition, convergent monotone methods generate solutions 
that automatically satisfy an entropy condition. This means that if the scheme is consistent and con- 
servative, then it not only generates a weak solution (Lax-Wendroff theorem [30]) but also a unique 
solution. Consequently, both smooth and discontinuous flow solutions are necessarily the physically 
correct solutions. We should also point out that these monotonicity-preserving schemes impose the 
total variation bounded (TVB) condition, which is a weaker condition than TVD. This is because 
TVD schemes are at most first-order accurate (see theorems of Godunov [31], Harten [32]). 

There are several approaches to building schemes that preserve monotonicity. One is to start with 
a low-order monotonic scheme and introduce anti-diffusion that is limited to ensure local monotonic- 
ity (e.g., FCT scheme). Another approach is to use a monotone interpolation method (e.g., MUSCL) 
to determine the primitive variables for reconstruction of the fluxes at the cell interfaces. Limiting is 
applied to the slopes used in the interpolation to impose monotonicity. A third way is to blend low- 
order and high-order diffusive terms and choose appropriate coefficients to ensure that the scheme 
is positive; and thus, the scheme is local extremum diminishing (LED). The concept of LED was 
introduced by Jameson [33,34] to provide a theoretical foundation for extending the non-oscillatory 
design principle (i.e., principle of nonincreasing maxima and nondecreasing minima) to multiple di- 
mensions, since the TVD property is a one-dimensional concept. In all of these methods the limiter 
must satisfy certain properties. In general, these are natural averaging properties and the property 
to ensure satisfaction of local monotonicity (e.g., LED). These schemes are also designed to pre- 
serve positivity. In addition, by relaxing the LED property (i.e., essentially LED), these schemes 
can be made higher order even at smooth extrema. Although most of the theory for non-oscillatory 
algorithms has been developed using scalar equations and constant coefficient one-dimensional sys- 
tems, these nonlinear high-resolution schemes have proved to be effective for the multidimensional 
nonlinear systems of fluid dynamics. 

So far we have discussed primarily second-order non-oscillatory methods, with some limited 
reference to the third-order scheme called PPM, which is based on a piecewise quadratic representa- 
tion of the solution. In general, higher order methods (order greater than two) rely upon techniques 
such as filtering and/or nonlinear coefficients (e.g., essentially non-oscillatory schemes) to prevent 
oscillations at extrema and to ensure stability (i.e., to prevent unbounded energy growth due to high- 
frequency error components). For example, Yee and Sjogreen [35] developed a class of high-order 
central difference schemes that apply a multistep filter method that includes high-order linear and 
nonlinear filters. The linear filter consists of a product of a sensor and either a spectral-like filter or 
a high-order centered dissipative operator. The nonlinear filter consists of a product of a different 
sensor and either an artificial compression method indicator or a wavelet sensor and the dissipative 
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part of a high-resolution shock-capturing algorithm. Another high-order approach for preventing 
spurious oscillations at extrema includes the essentially non-oscillatory (ENO) and weighted es- 
sentially non-oscillatory (WENO) schemes (see Shu [36]). The ENO schemes are Godunov-type 
schemes that achieve high-order spatial accuracy in smooth regions by a piecewise polynomial ap- 
proximation of the solution based on cell-averaged values. A nonlinear adaptive procedure is used 
to automatically choose the smoothest local approximation stencil, and thus, avoid the Gibbs phe- 
nomenon at discontinuities (i.e., generation of 0(1) spurious oscillations). Both the reconstruction 
and the resulting scheme are essentially non-oscillatory (i.e., almost TVD). The ENO schemes are 
essentially non-oscillatory because they allow small spurious oscillations that have size on the order 
of h r , where h is the mesh spacing and r is the approximation order of the scheme. Instead of using 
the optimal stencil, the WENO schemes use a convex combination of the candidate stencils with 
nonlinear adaptive coefficients (weights). In practice, both the ENO and WENO schemes generally 
produce solutions without oscillations. 

There are also hybrid high-order algorithms that use a WENO algorithm in the neighborhood of 
discontinuities and non-dissipative algorithms such as a tuned center difference method or spectral 
method (exponential convergence method) for smooth regions (e.g.. Hill and Pullin [37], Costa and 
Don et al. [38]). Another type of high-order scheme that is frequently used is the compact scheme, 
which is based upon a Pade differencing operator. Examples of such schemes are the centered 
sixth-order and fourth-order schemes using five and three points, respectively, employed by Visbal 
and Rizzeta [39] and Visbal and Gaitonde [40], They use a high-order low-pass spatial filtering 
technique to provide dissipation at the high modified wavenumbers (preventing energy build up at 
the high frequencies) when there are significant dispersion errors due to the spatial discretization. 
At discontinuities they use a Roe scheme with the MUSCL approach for reconstruction, which 
produces a third-order upwind-biased scheme in the case of uniform grids. Various types of detection 
mechanisms have been considered with these compact schemes, and they include: (1) Jameson-type 
pressure switch, (2) the artificial compression method, (3) multiresolution wavelets, and (4) WENO- 
type smoothness criteria. 

At this point we have simply identified different types of ILES schemes. Now, we need to be 
more specific about which ILES schemes are appropriate for indirectly introducing a subgrid scale 
model. Grinstein and Fureby [8] have demonstrated that the leading order truncation error terms for 
a certain class of ILES methods (i.e., methods having a flux function determined by blending of low- 
order and high-order fluxes with an appropriate limiter function) provide implicit SGS models that 
are similar in form to conventional (explicit) mixed SGS models. While a general theory for the ILES 
approach is yet to be established, the results obtained for fundamental turbulence problems such as 
isotropic turbulence decay and Taylor-Green vortex problems provide some support to the validity 
of the approach. Nevertheless, general validity can only be established by directly connecting the 
ILES methodology to the physical behavior of small scale turbulence. 

For both conventional LES (with explicit SGS models) and ILES (for which the SGS model is 
“implicitly” included through the numerical scheme), the magnitude of the numerical errors relative 
to the magnitude of the subgrid terms is critical in determining the accuracy of LES computations. 
There are two fundamental types of numerical errors: discretization (truncation) errors and alias- 
ing errors. The truncation errors are determined by the order of approximation of the terms in the 
governing flow equations. Aliasing errors are caused by the interaction of nonlinear terms, result- 
ing in high-frequency errors reappearing as low-frequency errors. For central differencing and pure 
upwind differencing the leading truncation error terms are dispersive and dissipative, respectively. 
Upwind biased schemes have both types of errors. In order to accurately represent the turbulent ki- 
netic energy spectrum the difference approximation of convective terms must have low dissipation. 
At the same time, the dissipation should be sufficiently large to suppress oscillations and ensure 
numerical stability. In addition, there is a need to have appropriate dissipation for shock capturing. 
There is also a need to have low dispersion errors. For any scheme being used for LES, these nu- 
merical errors must be monitored and controlled to ensure that they are sufficiently small so that the 
dominant errors in the simulation are due to SGS modeling. Only then will it be possible to effec- 
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tively improve the models and obtain a simulation capability for a broad range of flow conditions. 
Since many existing schemes use upwind biased spatial differencing for shock capturing, there is 
clearly a need to investigate and validate these schemes for LES computations. 

1.4 Related Work and Current Contribution 

Gamier et al. [41] investigated six shock capturing schemes for LES computations. These methods 
included the following: Jameson scheme with scalar dissipation. Roe scheme with the MUSCL ex- 
trapolation procedure and two minrnod limiter variations, and three different types of ENO schemes. 
Their spatial accuracy range was between second and fifth order, and either second-order or third- 
order TVD Runge-Kutta (RK) schemes were used for time integration. Each scheme was applied to 
decaying isotropic turbulence problems ignoring viscous effects, and Mach numbers were between 
0.2 and 1.0. Simulations were performed for mesh resolutions of 64 3 and 128 3 . The effects of both 
implicit and explicit SGS models were discussed. Gamier et al. [41] considered a scheme suitable 
for LES if it satisfied one of the following two conditions: (1) numerical dissipation is much smaller 
than the physical SGS dissipation, or (2) numerical dissipation mimics that of a SGS model. They 
concluded that none of the schemes tested satisfied either of these two conditions. 

Another study concerning the behavior of high-resolution and high-order methods for implicit 
LES was performed by Thornber et al. [42], Their investigation, at low Mach number, was also done 
for homogeneous decaying turbulence. Computing on mesh densities between 32 3 and 256 3 , they 
considered seven different schemes with spatial accuracy ranging from second to ninth order. In 
all simulations the Euler equations were solved with finite-volume methods. For higher order spa- 
tial accuracy, cell interface variables were determined by either the MUSCL or WENO approach. 
Third-order RK schemes were used to advance the solution in time. With grids exceeding 32 3 or 
with numerical methods having spatial accuracy higher than third order, their results indicated that 
the growth of large scales and the dissipation of kinetic energy could be captured. Based on en- 
ergy spectra, Thornber et al. also concluded that all the upwind-biased schemes evaluated were too 
dissipative at the high wavenumbers, resulting in an early onset of decay in the inertial range. We 
should point out that with this approach the question does arise concerning an appropriate repre- 
sentation of an inertial subrange in the absence of physical viscous effects (i.e., the physical effects 
when Re = oo are clearly not the same as for high Re). Thornber et al. also observed that due to 
differences in the effective cut-off wavenumber, methods with fifth-order or higher spatial accuracy 
performed better on a 128 3 grid than those with second-order accuracy on a 256 3 grid. 

The primary objective of this report is to compare the behavior of representative second-order 
finite-volume schemes and high-order finite-difference WENO-based schemes when applied to the 
problem of isotropic decay of turbulence. The second-order schemes have either central or third- 
order upwind biased spatial differencing, and for each scheme second-order temporal differenc- 
ing is used. Both explicit and implicit LES computations are performed with these methods. For 
the WENO schemes we investigate fifth-order, seventh-order, and ninth-order spatial accuracy for 
ILES. A RK scheme, which is third-order TVD, is used for time advancement. There is a recogni- 
tion (Martin et al. [43], Taylor and Martin [44]) that the standard fifth-order WENO scheme with 
Lax-Friedrichs flux splitting is too dissipative at the higher wavenumbers, as clearly revealed in the 
energy spectrum. Here, the effect of increasing the order of approximation is discussed. In the simu- 
lations, the filtered Navier-Stokes equations are solved. Since the Kolmogorov scale is significantly 
smaller than the mesh size. Gamier et al. [41] and Thornber et al. [42] indicated either indirectly 
or directly that the viscous effects are negligible. However, in any practical simulation the viscous 
terms and their effects are present, regardless of the spatial resolution. This certainly suggests that 
including these terms is necessary to obtain an effective measure of the wave cut-off number in the 
energy spectrum. Moreover, as discussed by Ishihara [45], the turbulence statistics have significant 
Re dependencies, even at high values of Re. 

In the initial part of the paper the governing equations are presented. Then the numerical methods 
are described. This is followed by the results section, which begins by describing the formulation of 


7 



the initial solution and providing definitions of the turbulence statistics. In the next part the results 
from LES simulations are presented. The turbulence results include kinetic energy, dissipation of 
energy, and skewness. Vorticity and Q-invariant contours are also shown. The initial peak of the 
energy dissipation during the energy redistribution process clearly reveals the improved accuracy 
achieved with the higher order WENO schemes. The computed energy spectrum is used to determine 
the dissipative behavior of the numerical methods over the frequency range. Results with the various 
numerical schemes are compared with each other. In addition, there are some comparisons with 
results from a spectral scheme. 


2 Governing Equations 

2.1 Direct Numerical Simulation (DNS) 


Written in strong conservation form, the compressible form of the Navier-Stokes equations is: 

dp_ dpu k = 
dt dxk 


(1) 


dpui dpuiUk 


dt. 


dxk 


dp 

dxi 


d 

dx k 


dui duk 2 duj 

dxk dxi d dXj 


( 2 ) 


dpE d[(pE + p)u k ] _ d 
dt dx k dx k 


j dT 

dx k 




dui 

dx k 


du k 2 duj 

~ 7,0ki 


dxt 


dxn 


( 3 ) 


where the coefficient of thermal conductivity k! is taken to be k' = c p p/ Pr, and the variable E 
represents the specific total energy E = e + ( u 2 + v 2 + w 2 )/2. For a perfect gas, e = c v T, 
c v = R /( 7 — 1), 7 = c p /c v , and the equation of state is p = pRT , or: 


P=( 7 - 1 ) 


pE — i p ( u 2 + v 2 + w 2 ) 


( 4 ) 


By definition, DNS is a numerical computation of the Navier-Stokes equations that resolves all 
relevant spatial and temporal scales in a flow field. This means that the grid needs to be fine enough 
to resolve features on the order of the Kolmogorov dissipation length scale: 

p=(v 3 /e) z (5) 

where v is the kinematic viscosity and e is the dissipation rate. Similarly, the time step needs to be 
fine enough to resolve temporal dynamics on the order of the Kolmogorov time scale: 

T=(v/e)f ( 6 ) 


It is impossible today to achieve these resolutions at reasonably high Reynolds numbers. Successful 
low Reynolds number simulations require very low-dissipation schemes to resolve the smallest- 
scale features with good fidelity. The spatial order of accuracy of the scheme also influences grid 
requirements. Additional dissipation is required to stabilize some numerical schemes; for example, 
compact schemes are often stabilized through the use of additional filtering [46]. Note that an under- 
resolved DNS simulation can exhibit time-dependent features that appear to the eye like an adequate 
DNS simulation, but if the smallest space and time scales are not resolved, some flow physics are 
being omitted and its validity remains unclear. 



2.2 Large Eddy Simulation (LES) 

In LES, the Navier-Stokes equations are filtered via a low-pass filter defined as a convolution prod- 
uct: 


</>(£, 




<K£t')G(x 




t')dt'd 3 £ 


(7) 


where < p(x,t ) is the resolved part of a space-time variable <f>(x,t). For the compressible Navier- 
Stokes equations, the equations contain both ordinary and Favre-filtered variables (see, e.g., Knight 
et al. [47]). The Favre filter is defined by: 


/ = — ( 8 ) 

P 

where the overbar denotes the ordinary spatially filtered variable, and the tilde denotes the Favre- 
filtered variable. After filtering, the compressible governing equations (Eqs. (1) - (3)) become: 


dpu k 
dt dxk 


(9) 


dpiii dptiiUk dp dT lk 


dt 


dxi 


dxj dxi 


dpE d[(pE + p)u k ] dHk 


dt 


dxk 


dxk 


( 10 ) 


( 11 ) 


where the total stress tensor is 7 = J>Ti k — o’ lk , and the sum of the heat flux plus work done 
by the stresses is Tt k = Qk + n'dT/dxk + EikEi with k' typically expressed as c p p,/ Pr, and a lk 
approximated as: 


/ duj 
\dx k 


du k 

dxi 



( 12 ) 


As a result of the filtering, the effect of the filtered stress kinetic energy (k = th/2) should be 
included in the definition of the total energy (Knight et al. [47]), that is, J>E = pe + p(u 2 + v 2 + 
w 2 )/2 + pk. Then, the equation of state becomes: 


P=( 7-1) 




v 2 + w 2 ) — pk 


(13) 


However, for lower speed flows, the effect of k is usually negligible, and many researchers ignore it 
in the definition of energy and equation of state. 

As a result of the filtering, two unknown quantities emerge: the subgrid scale stress T,k and the 
heat flux Q k : 


Tik = UiU k - UiUk (14) 

Qk = -c p p (u k T - u k Tj (15) 

The quantities in Eqs. (14) and (15) must be modeled. As discussed in Sagaut [9], there are two 
basic methods: explicit modeling and implicit modeling. In explicit modeling, a subgrid model 
is explicitly introduced. In implicit modeling (ILES), no additional eddy viscosity is introduced 
into the governing equations, but the numerical method is selected such that the numerical error 
fulfills desired properties and effectively acts like a subgrid model. With no explicit model present, 
ILES appears to be functionally similar (from an implementation point of view) to under-resolved 
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DNS; neither method resolves the finest scales, and both require enough inherent dissipation in the 
numerical scheme to prevent the non-physical build-up of energy at the smallest resolved scales. 

The commonly used Smagorinsky subgrid scale model [2] is an explicit SGS model, which 
utilizes the the Boussinesq assumption for the SGS stress tensor: 




&ik 0 


2 f , 
+ gOjfeK 


G6) 


where the v t in the Smagorinsky model is given by: 
v t = (c s A) 2 |S| 


(17) 


with |S| = yj 2 SijSij, Sij = (dv,i/dxj + duj /dxi) /2, and A is typically defined by some measure 

of the local grid spacing, such as A = (AxAyAz)^ 1 ^ 3 ^ . The c s term is a constant for the original 
Smagorinsky model, commonly set to approximately 0.1. 

Note that in most subgrid scale models, neither the formal filter function nor the filter width are 
explicitly defined. Instead, the models have a built-in filter, related to the local grid spacing: for the 
Smagorinsky model it is c s A. This built-in filter controls the size of the smallest locally resolved 
flow structures. It can be difficult to choose an optimum built-in filter width, which has been noted to 
be flow dependent. This difficulty has been one reason for the development and subsequent success 
of the widely used dynamic model of Germano et al. [3], This model dynamically computes a 
variable c s term, rather than setting it to a constant value. 

For compressible flows, an SGS viscosity model for the subgrid-scale heat flux is (see Urbin and 
Knight [48]): 


_ pc p vt dT 
Pr t dx k 


(18) 


3 Numerical Methods 

For the current problem of decay of isotropic turbulence, two compressible CFD codes were em- 
ployed, the NASA production computer code CFL3D [49] and an in-house NASA research code that 
employs the weighted essentially non-oscillatory (WENO) scheme [36]. Some of the characteristics 
of these numerical methods will be described next. 


3.1 Upwinding and Central Differencing in CFL3D 

The CFL3D code is a structured-grid upwind multi-zone CFD code that solves the generalized thin- 
layer or full Navier-Stokes equations [49], The code can employ local time-step scaling, grid se- 
quencing and multigrid to accelerate convergence to steady state. A second-order time-accurate 
mode (using standard backward differentiation formula (BDF)) is used in this study with pseudo- 
time stepping plus multigrid, and an implicit approximate factorization method to advance the solu- 
tion to the next physical time step. For the current study, a 3-level W-cycle multigrid strategy was 
employed, along with 30 sub-iterations (pseudo-time steps) per time step. This number of subitera- 
tions was enough to drive the L 2 -norm of the subiteration density residual down approximately 14 
orders of magnitude, to near machine zero for the current problem. 

The CFL3D code employs a cell-centered finite-volume method. It uses the so-called “kappa 
scheme” in a MUSCL approach [50] for its spatial differencing for the convective and pressure 
terms, and second-order central differencing for the viscous terms. Roe’s flux difference-splitting 
(FDS) method [51] is used to obtain fluxes at the cell faces. In the kappa scheme, the reconstruction 
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at cell interfaces is accomplished via: 


U L = U 1 +^[(1-k)(U 1 -U 1 . 1 ) + (1 + k)(U 1+1 -U 1 )} (19) 

U R = U i+1 -^[(l + K)(Ui +1 -Ui) + (l-K)(U i+ 2-U i+ 1)] (20) 

where Ui is the state vector at cell center i, and the subscripts L and R refer to the left and right 
states for a local approximate Riemann solver. The “standard” scheme for CFL3D in RANS appli- 
cations uses k = 1/3, which is a third-order upwind-biased stencil on a uniformly spaced Cartesian 
grid [52]. For the LES applications in the current paper, upwinding was believed to be too dissipa- 
tive, so for most of the applications the central-differencing option k = 1 was used. When k = 1, 

U L = U R =^(Ui + U i+ 1) (21) 

The Roe-type numerical flux is written as 

fi+ 1/2 = + /«) “ 2 l^li+1/2 ( Ur ~ Ul ) ( 22 ) 

Since Ur = Ur , the dissipative term in Eq. (22) vanishes, and the flux balance for a given cell 
yields a standard central difference approximation with second-order spatial accuracy. In general, 
central differencing with no additional dissipation can yield numerical problems such as odd-even 
point decoupling. However, for the particular case of isotropic turbulence decay considered here, 
the solutions on the grid sizes considered were stabilized sufficiently by the viscosity of the fluid 
(although not shown, computations attempted on a very fine grid with 256 3 cells did yield odd-even 
point decoupling behavior). 

The standard Smagorinsky model has been coded into CFL3D, with the constant coefficient c s 
specified by the user. When c s is set to 0, the method becomes an ILES type, with no additional 
subgrid eddy viscosity model added. 

3.2 Weighted Essentially Non-Oscillatory (WENO) Scheme 

In this section, for completeness, we present the elements of the weighted essentially non-oscillatory 
(WENO) scheme and the implementation of this scheme in the present work. To simplify the de- 
scription of the scheme, we consider a one-dimensional scalar conservation equation 

du df{u) _ df v (u x ) 

dt dx dx ( ’ 

where u — u(x, t) is a conserved quantity, f(u) and f v (u x ) are inviscid and viscous flux functions 
with u x denoting the derivative du/dx, and the independent variables x and t represent space and 
time, respectively. By taking this approach there is no compromise of the generality of the discus- 
sion since the same spatial approximation techniques can be applied in each coordinate direction of 
a multidimensional fluid dynamics problem. We obtain a semi-discrete form of Eq. (23) by consid- 
ering a one-dimensional domain with uniform spacing Ax. The unknowns are located at the points 
Xi = iAx, i = 0, . . . , N, and cell boundaries are at x i+1 /2 = x* + Ax/2. As seen in Fig. 1, the 
grid points are also cell centers. The semi-discrete form of Eq. (23) is given by 

duj(t) = _ df_ + dp_ 
dt dx x=x , dx 

where -«,;(£) is a numerical approximation of the point value u(xi,t ). To construct a conservative 
difference scheme for this system of ordinary differential equations, the fluxes must be evaluated 
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at the interfaces of the cells to the desired order of accuracy. Since the fifth-order WENO scheme 
is frequently used in applications, we present the details of the discretization for the derivative of 
/(it) using a fifth-order finite-difference form of this scheme, following the work of Shu [36], Jiang 
and Shu [53], Henrick et al. [54], and Borges et al. [55], After discussing the approximation of 
convective term, we then consider different methods for discretizing the viscous (i.e., diffusion) 
term, giving the advantages and disadvantages of each method. Next, we provide a brief discussion 
of the extension to even higher order based on the paper by Balsara and Shu [56]. In the last part of 
this subsection we describe the time integration method. 

3.2.1 Approximating the Inviscid Term 

This section contain three parts. In the first part we consider the linear form of the discretization for 
the derivative of the convective (inviscid) flux function. This form is dominant in smooth regions 
of the domain, and it provides a natural foundation for the development of the nonlinear WENO 
scheme. The final part of this section discusses the flux splitting used to account for different signs 
of the flux Jacobian. 


3.2.1a Linear Form of 5th Order WENO Scheme 


Consider the inviscid flux function. One can satisfy the conservative property by defining a numeri- 
cal flux function h(x) implicitly as 



so that the spatial derivative of f(x) is exact; and thus. 


df 

dx 




+ 1/2 ~ hi- 1/2 


(26) 


where hi±\/ 2 = ^(**± 1 / 2 )- The numerical flux function h i±1 / 2 can be approximated by a flux 
function fi±\/ 2 that is defined by a polynomial of sufficient degree to obtain the desired order of 
accuracy. Then 


d / 

dx 


^:(/*+i /2 ~ /*- 1 / 2 ) 


and 


(27) 


fi+ 1/2 — M^+l/ 2 ) + 0(Ax r ) 


(28) 


where r is the desired order of accuracy. The numerical flux function fi+ 1/2 depends on the discrete 
values of the physical flux / at the points which determine the stencil of the approximation. 

For the WENO scheme a fifth-order accurate flux fi+ 1/2 is computed as a convex combination 
(i.e., weighted average) of possible third-order approximations. In Fig. 2 the three candidate stencils 
are shown. The weighted average of the third-order stencils can be expressed as 

r— 1 

hi+ 1/2 ~ fi+ 1/2 = Wfc /? i +i /2 (29) 

k - 0 

where r = 3 and +/, is the weight corresponding to stencil 5+. To obtain the third-order approxima- 
tions for each f k (x), we represent each approximate numerical flux function with a second-degree 
polynomial, and thus 

h(x) « f(x) = ao + a\x + a 2 x 2 (30) 
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with the undetermined coefficients a+, k = 
integrating we obtain 

f(x) = a 0 + a\X + a 2 



0, 1, 2. Substituting f(x) for h(x) in Eq. (25) and 


(31) 


The coefficients a*, of the polynomial are determined from values of f(x) at nodes of the grid. Con- 
sider the grid cell with boundaries at Xj_i/ 2 and x i+1 / 2 - To obtain the a*, so that the flux of Eq. (29) 
produces a fifth-order upstream-centered scheme at the interface x i+1 / 2 , the values of f(x) are 
required at the following three sets of nodes: {xj_ 2 , Xi-i, Xi}, Xi, x, + i}, {xi, tCj+i, Xj +2 }. 

Evaluating the polynomial of Eq. (31) for each set of nodes gives three algebraic equations that must 
be solved to obtain the coefficients of the interpolating polynomial. For the evaluations the value of 
Xi is zero. The numerical flux functions for the three sets of points are 


f°(x) 

;v) 

f\x) 


—fi-2 + 2/j-! + 23 fj /7i- 2 -4/i-! + 3/, ; 

24 H 


—fj-i + 26 fj + fj + 1 
24 

(23 fi + 2 f i+ i — fi+ 2 ) 
24 


2 Ax 
fi+i fi— 1 


fi-2 ~ 2/j_i + 


2Ax 2 

fi-l ~ 2fi + fi+l\ 2 


2Ax 

—3fi + 4/i + i — fi + 2 

2 Ax 


2Ax 2 


(32) 


fi ~ 2fi+l + /i+2 \ x 2 
2Ax 2 


At i + 1/2 (x i+ i /2 = Ax/2) the fluxes of Eq. (32) become 

fi+ 1/2 = g(2/i-2 — 7/t-i + H/i) 

fl+ 1/ 2 = ^(-fi-i+5fi + 2f i+1 ) (33) 

fi+ 1/2 = g (2/t + 5/ i+ i — fi+ 2 ) 

One can easily show (see Henrick et al. [54]) that the fifth-order flux at .x', + | / 2 for the upstream- 
centered scheme is given by 

fi + !/2 = ^ (Vi-2 - 13/j-i + 47/,; + 27 fi+i - 3 f i+2 ) (34) 

Remark: The third-order fluxes of Eq. (33) and the fifth-order flux of Eq. (34) can also be determined 
by using Table 2.1 in the 1997 paper by Shu [36]. 

By substituting the third-order approximations of Eq. (33) into Eq. (29) and equating corre- 
sponding terms, we obtain the weights 


cuo — 1/10, cui — 6/10, 1 x 2 — 3/10 


(35) 


which give the desired fifth-order flux function in regions where the solution is smooth. The weights 
sum to unity, and they are called the ideal weights. To obtain the flux fi- 1/2 for approximating the 
first derivative of the inviscid flux function /, the stencil for fi+ 1/2 is shifted to the left by one mesh 
point. Substituting for the fluxes in Eq. (27) gives 

« 2/i— 3 + 15 /i— 2-60/i — 1 + 20/i + 30/ m - 3/ i+2 ) (36) 

= f'(x) + 0(Ax 5 ) 


V 

dx 


Even though the fluxes at Xi _\/2 and x i+1 /2 are only fifth-order accurate and the flux balance of 
Eq. (27) is divided by Ax, the resulting flux derivative is fifth-order accurate because the two stencils 
approximating the fluxes only differ by a grid point shift, resulting in a canceling of the lowest order 
truncation error terms. 
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3.2.1b 5th Order Non-Oscillatory Weights 

So far we have only considered a linear scheme. The approximation of the inviscid flux derivative 
given by Eq. (36) requires a smooth solution for high-order accuracy. If a discontinuity such as a 
shock wave occurs in the flow field, oscillations will occur in the discrete solution. The full nonlinear 
WENO scheme can be designed to give comparable high-order accuracy in smooth regions to that 
of the linear scheme and prevent spurious oscillations at discontinuities. This is accomplished by 
introducing appropriate nonlinear weights in Eq. (29) such that the weight of each stencil depends 
upon its smoothness. In the nonlinear scheme, stencils that cross a shock wave are automatically 
assigned an essentially zero weight, producing a one-sided upwind difference approximation. Such 
an approach is similar to what is done with second-order schemes, which use a limiter for shock 
capturing that produces a first-order upwind difference at shock waves. One can view the nonlinear 
weights of the WENO scheme as playing the role of a limiter, and by their action, the scheme is 
essentially non-oscillatory. In the neighborhood of a discontinuity the order of accuracy is reduced. 
That is, if r is the order of accuracy of the weighted fluxes, then in regions where the solution is 
smooth the order of accuracy of the scheme is 2 r — 1 and in the vicinity of a shock wave the order 
of accuracy is reduced to r. 

Jiang and Shu [53] define the non-oscillatory weights as 


Wfc 





dk 

(/3k + s w ) p 


(37) 


where the are normalized weights, dk are the ideal weights, and the parameter s w is a small 
positive number used to bound the weights (i.e., to prevent division by zero). For r = 3 the exponent 
p = 2 is chosen, since it proves to be adequate to obtain essentially non-oscillatory approximations. 
The /3k is a smoothness measure of the flux function on the stencil ,S/ , and it is defined by 


r— 1 


f3k = ^2 


21—1 


l - 1 


r Xi + 1/2 ( d l f k 

'Xi- 1/2 


dx l 


dx 


(38) 


Since f k = oo + a\X + a 2 X 2 , the equation for f3k becomes 
(3 k = a 2 Ax 2 + ^a 2 Ax 4 

o 


(39) 


when the stencils are centered around Xi = 0. Substituting for the coefficients « | and 02 corre- 
sponding to each value of k we obtain the smoothness indicators 


Ao = Y^(fi-2 — 2/i-l + fi) 2 + -(/i— 2 — 4/,_! + 3/i)" 

At = Y^(fi-l ~ 2fi + fi+l) 2 + ^(fi+1 - fi-l) 2 

A 2 = ^(/i — 2/i+i + fi+ 2) 2 + -(3/j — 4/i+i + fi+ 2) 2 


(40) 


which have Taylor series expansions at Xi given by 


/ 3 0 = /'Ax 2 + (g/f - |///r) ax 4 - (y /r/r - a * 5 + o(Ax 6 ) 

At = /'Ax 2 + + Ax 4 + 0(Ax 6 ) 


( 41 ) 
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A 2 = /'Aa- 2 + -/r - -flfl' Ax 4 + -/"/'" - -fir Ax 5 + 0(Ax b ) 


13 
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with the prime denoting differentiation (e.g., f" = d 3 f/dx 3 ). Appropriate smoothness indicators 
have the following behavior. In regions where the solution is smooth, the 0 k are small and about the 
same size, resulting in weights u>k that are close to the ideal weights d k - If the approximation stencil 
contains a discontinuity, then the associated smoothness measure 0 k is 0(1), and the corresponding 
weight LOk is small relative to the other weights. 

When the solution is smooth, a requirement in choosing the appropriate 0k’ s is that the corre- 
sponding Wfc’s approach the ideal weights at a sufficiently fast convergence rate. Henrick et al. [54] 
and Borges et al. [55] derive the necessary and sufficient conditions to ensure the design order of 
accuracy in smooth regions. These conditions are as follows: 

r— 1 

- d k) = 0(Ax 2r ) (42) 

1=1 

r—1 

E A k{^t - ) = O(A^) (43) 

i=i 

u±-d k = 0(Ax r - 1 ) (44) 

where .4 /, is the coefficient (independent of Ax) of the leading term in a Taylor series expansion 
fOT fi± 1/2. and the superscript ± refers to this flux function. As indicated by Borges et al. the first 

condition is always satisfied, since ^ = Sfc=o^ fe - A sufficient condition for fifth-order 

convergence is given by 

u k - d k = 0(Ax 3 ) (45) 


Following the analysis of Henrick et al. for the smoothness indicator defined by Eq. (40), one can 
show that if 0 k can be expressed as 

0k = D[1 + 0(Ax r_1 )] (46) 

where D is a non-zero constant independent of k, then a u k = d k + 0(Ax r ). When r = 3 we 
readily see from the Taylor series expansions for the 0 k s given in Eq. (41) that Eq. (46) is satisfied, 
and thus, to k = d k + 0(Ax 2 ). Borges et al. emphasize that the condition of Eq. (43) must be 
satisfied to achieve the design order of convergence. 

At critical points, which are points where the first or higher derivatives of the function / vanish, 
these conditions are not necessarily satisfied using the 0 k given in Eq. (40). For example, if /' = 0, 
we find from the Taylor series expansions of Eq. (41) that 0k = D[ 1 + O(Ax)] for k = 0, 2, and 
thus, u>k = d k 4-0 ( Ax j, resulting in third-order convergence. In order to overcome this deterioration 
in the convergence of the scheme due to critical points, Borges et al. introduced a different definition 
of the measures of smoothness 0 k given by 


0 


Z 

k 


f 0k T &W \ 

\ 0 k + T 5 + £ w ) 


k = 0, 1, 2 


(47) 


where 75 = \ 0 O — 0 2 1, and the superscript 2 indicates that the quantity, a smoothness measure in this 
case, is associated with the WENO-Z scheme of Borges et al. The corresponding weights u>% are 
defined as 


, .* _ a l n z _ dk _ , 

w fc \r— 1 z ’ a k o z dk 

J2k = 0 a l Pk 


1 + 


T~5 


0 k + 


k = 0, 1, 2 


(48) 


The exponent q > 1 increases the rate of convergence to the ideal weights for smooth regions of the 
solution, and it accelerates the convergence rate to zero in the vicinity of discontinuities. 

In practice, the question arises as to whether or not the uj k are sufficiently good approximations 
of the ideal weights for smooth parts of the solution. The appropriate scaling of the smoothness 
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indicators /3k is essential to ensure that they do not dominate the parameter s w in smooth regions and 
produce weights that will result in too much numerical dissipation. This suggests that a given scaling 
of [3k can have an effect on the proper choice of e w . At the same time, caution must be exercised to 
make certain that the choice of e w does not compromise the non-oscillatory behavior of the WENO 
scheme. In numerical computations we have observed that if the parameter e w is too small relative 
to the smoothness measures defined by either Eq. (40) or Eq. (47), the WENO scheme will have too 
much numerical dissipation. Other difficulties can also occur when e w is not sufficiently large. Shen 
and Zha [57] found that when e w is too small (e.g., 10~ b ) relative to /3k, convergence and stability 
problems can occur due to small disturbances in the [3k . These disturbances result in oscillations of 
the nonlinear weights which prevent convergence to the ideal (optimal) weights in smooth regions. 
Thus, there are increases in the numerical dissipation. Shen and Zha considered a seventh-order 
accurate scheme with the smoothness indicators of Borges et al. [55] and solved either the two- 
dimensional Euler or Navier-Stokes equations for subsonic and transonic steady flows. By setting 
e w = 10 -2 , Shen and Zha prevented convergence stall and showed improved accuracy. 

Oscillations in the weights are not only caused by flow discontinuities. As pointed out by Shen 
and Zha [58], strong variations in the smoothness measures (3k can occur due to nonuniformities in 
the computational mesh. If e w is very small (e.g., 10” 6 ), the (3k will not always be less than e w , 
even in smooth regions of the flow field. Here again, by increasing the magnitude of e w one can 
avoid oscillations in the nonlinear weights. 

In the present work we apply the weights of Eq. (47). For comparison purposes we also consider 
the weights of Eq. (37). The influence of the parameter e w on accuracy is shown in the results 
section. 


3.2.1c Flux Splitting 

Throughout the previous discussion we have assumed that the flux Jacobian (df /du) is positive. 
In general, there is a need to allow for a sign change in the flux Jacobian. This is usually done by 
splitting the flux into positive and negative parts. In the present formulation for the WENO scheme 
a global Lax-Friedrichs splitting is used. Thus, 

f ± W) = ^(f(u)±au) (49) 

where a = max\df /du\, and the maximum is taken over the complete solution field. A local 
splitting has also been considered. For the decay of turbulence problem considered herein, the local 
splitting did not have any significant effect on the results. 


3.2.2 Approximating the Viscous Term 

Here we discuss the approximation of the physical diffusion terms in the Navier-Stokes equations 
by considering the derivative of the viscous flux function f v {u x ) in Eq. (23). The derivative of the 
viscous flux function is defined by 


df v (x ) d [ , . du(x ) 


(50) 


where in general, the coefficient a( x) = a[u(x),x\. Since we want to construct a conservative 
difference scheme, this derivative is approximated by 

df V ( X ) ~ fl’ +1/2 - fj- 1/2 
dx Ax 

Thus, one needs to determine the interface approximations of a(x) and u(x) at Xj+i /2 and Xi_i/ 2 - 
With the fifth-order WENO scheme, the diffusion term is approximated with a conservative fourth- 
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order finite difference. When evaluating possible ways to determine these discrete interface func- 
tions, consideration must be given to not only accuracy and stability but also damping of all error 
modes. We now consider two ways to compute a^± 1/2 and u iz j =1 / 2 - 

One possible method of discretizing the diffusion term is to first compute the discrete f v (x) at 
the solution points. A fourth-order approximation of f v (x ) at x t , an interior point of the domain, is 
given by 

U v )i = 12 Ax - + 8ui+1 ~ u *+ 2 ] (52) 

Then the approximation for the derivative of f v (x) is evaluated when constructing the approximation 
for the derivative of f(x). For the fifth-order linear form of the WENO scheme (i.e., ideal weights) 
we obtain 

[-2(/")i-3 + 15(r )i -2 - 60(n,_ 1 + 20 (/+ (53) 

+ 30(r) i+1 -3CT) i+2 ] 

— — — — 2 ^ a i-3 u i-5 + (16a^_3 + l5a,i-2) u i-4 

72(JAx 

— (120ai_2 + 60ai_i)'Ui_ 3 + (— 16ai_ 3 + 480a^_i + 20 ai)ui-2 

H- (2fl^— 3 H - 120a^_ 2 — 160a^ H- 30o.2+i)^i— 1 

— ( \hcii—2 + 480ai_i + 240a^_(_i + 3cii+2) u i 
+ (60(3.2—1 + 160(32 + 24^+2)332+1 + ( — 20(32 + 240(32+1)332+2 

— (30a2+i + 24 ^+ 2 ) 332+3 + 3 ^+ 2332 + 4 ] 

If a(x) and u(x) are smooth functions, then the approximation of Eq. (53) is fourth-order accurate. 
While this method for differencing the derivative is convenient for implementation and operation 
count, it does have the undesirable consequence that the 7r-mode is not damped. One can easily 
show this by treating the discrete product of a(x) and u(x) as a single variable, Fourier transforming 
the derivative, and evaluating the transformed discrete operator at a Fourier angle equal to ir. 

An alternative method that provides a compact difference stencil and ensures damping of the 
7r-mode is given by Kamakoti and Pantano [59]. This approximation enforces discrete conservation. 
Let B be a matrix of stencil coefficients for approximating the viscous term. Then, the flux form is 
given by 


dr 

dx 


df v (x) 


dx 


Ax' 


P P fV fV 

1 D _ 32+1/2 Ji- 1/2 

2 / J / J 3l>772+p+l ? n+p+l ^3+m332+n — 


m=—p n=—p 


Ax 


(54) 


where 


fi+ 1/2 ^ ^ ' 'y ^dm+p,n+p a i+rri u i+n 


(55) 


m=— p+1 n=— p+1 


with M denoting the matrix of the unknown coefficients for the viscous flux to construct the differ- 
ence stencil. The widths of the flux stencil and the stencil for the variable coefficient second-order 
derivative are 2 p and 2 p + 1, respectively. One uses Taylor series expansions to determine the 
elements of the matrix M. 

If we expand Eq. (55) for a fourth-order approximation, we obtain after rearranging terms 

fi+ 1/2 = ^[(MnOi— i + Mi\ai + M^iai+i + + (56) 

(Mi2a,_i + M22CL1 + M^di+l + M^20’i+‘2)Ui + 

(Mi 30 j_i + M 2 30 j + M 33 aj_|_ 1 + M/^ai + 2)ui + \ + 

+ M24di + + M44ai+2)'Ui+2] 
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The stencil matrix M is given by 


1 

8 


1 J_ 
6 24 


M = 



i I _I 

24 6 8 J 


( 57 ) 


The resulting approximation of Eq. (54) for the viscous term is compact, requiring just five points 
for fourth-order accuracy. Thus, it has a reduced error relative to the approximation of Eq. (53), and 
it has the advantage that the 7r-mode is damped. 

The discretization method used to obtain Eq. (53) was employed in the original ENO (Essentially 
Non-oscillatory) code of Atkins [60] that was modified by Balakumar et al. [61] to create a WENO 
version. While the 7r-mode will be damped when there is deviation between the WENO weights and 
the ideal weights, this mode may not be adequately damped in regions where the diffusion terms 
dominate or the ideal weights are recovered. For this reason, along with reducing the truncation 
error and ensuring that this error is purely dissipative, the discretization method for Eq. (53) should 
be replaced by that for Eq. (54). 

In this initial work we have retained the discretization method for the diffusion terms used by 
Atkins [60] and Balakumar et al. [61]. For the WENO scheme the polynomial that interpolates 
f v (a:) at the solution points to the cell interfaces requires the smoothness indicators of Eq. (40) to 
determine the appropriate non-oscillatory stencil weights. We should mention that periodic bound- 
ary conditions are imposed for the isotropic turbulence problem considered in this paper. This is 
accomplished in the WENO code by introducing additional points outside the physical domain. 


3.2.3 Higher Order Schemes 

Higher order WENO schemes can be derived by following the steps presented in the first two subsec- 
tions of Section 3.2.1 for the fifth-order scheme, starting with Eq. (29) and the appropriate r. For the 
seventh-order and ninth-order schemes the fluxes in the summation are fourth-order and fifth-order, 
respectively. Thus, each flux f^ +1 / 9 for the seventh-order and ninth-order schemes is determined 
by third-degree and fourth-degree polynomials, respectively. The stencils for theses schemes that 
define the corresponding WENO flux f i+ 1 / 2 , the ideal weights c//, ; , and the smoothness indicators 
Pk are presented in the paper by Balsara and Shu [56]. For completeness, these quantities are in- 
cluded in Appendix A. It should be noted that the smoothness indicators in Appendix A include the 
scaling that arises when using Eq. (38), whereas those of Balsara and Shu do not. For extremely 
small values of the e w in Eq. (37) and Eq. (48) this does not matter. However, it does have an effect 
for larger values of e w such as 10~ 2 . 

For the seventh-order and ninth-order schemes the /”(*) of the diffusion term is discretized at 
interior points with sixth-order and eighth-order central differences, respectively. Stencils for these 
approximations can be obtained using Table 2.1 of Shu [36], which is supplemented by Table Al in 
Appendix A to include eighth-order. The sixth-order approximation is given by 

( f V )i = QQAx ^- Ui ~ 3 + 9u?_2 “ 45ui -! + 45w i+ 1 _ 9u i+2 + Ui+ 3 ) (58) 

and the eighth-order one is given by 

(Hi = (3ui_4 - 32uj_ 3 + 168u,;_2 - 672^ + 672 u i+1 (59) 

- 168iq + 2 + 32w i+3 - 3it i+4 ) 

The interpolating polynomials for obtaining the viscous fluxes at the cell interfaces are constructed 
with the fluxes and WENO weights given in Appendix A. 
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3.2.4 Time-Stepping Scheme 

To advance the solution in time, a third-order Runge-Kutta (RK) scheme that satisfies the total 
variation diminishing (TVD) property is used. A scheme is said to be TVD if 


TV(u n+1 ) < TV(u n ), TV = 

E K+l - U i\ 

i 

(60) 

The RK scheme has three stages given by 


u « = 

u n + A tCu n , 


u & = 

j(3u" + u (1) + 

(61) 

u = 

\(u n + 2u^ +2 AtCu^), 

O 


u n+1 = 

U ( 3 ) 

(62) 


where C is the discrete spatial operator, and the superscript n denotes time level (time t = nAt). 

4 Numerical Results 

4.1 Problem Definition 

4.1.1 Decay of Isotropic Turbulence (DOIT) 

In the classical problem of decay of isotropic turbulence, an initial turbulent field is allowed to 
decay in a “box-like” domain with periodicity imposed on all six of its sides. The turbulent statistics 
are collected and compared with experiments or with other numerical studies. Particular focus is 
often placed on the decay of the total turbulent kinetic energy with time, the initial growth followed 
by subsequent decay of the total turbulent dissipation with time, and the energy spectrum at given 
instants in time. 

This classical problem is a good discriminator for SGS models and numerical schemes: poor 
models or overly-dissipative schemes tend to diffuse the energy-containing structures, and conse- 
quently miss many of the characteristic behaviors that can be accurately captured by well-resolved 
DNS. See, for example, the experimental results of Compte-Bellot and Corrsin [62], and the nu- 
merical studies of Clark et al. [63], Vreman et al. [64], Mansour and Wray [65], Spyropoulos and 
Blaisdell [66], Seror et al. [67], Hughes et al. [68], Hansen et al. [69], and Thornber et al. [42]. 

In the current study, the domain uses side lengths of 27r. Various Cartesian grid sizes are em- 
ployed, ranging from the coarsest with 33 s grid points (32 s grid cells) to the finest with 193 s grid 
points (192 3 grid cells). 

To evaluate the numerical solutions to this problem, we use DNS results computed with a spec- 
tral scheme. Spectral solutions provide a good measure of accuracy due to their low discretization 
errors. They have the additional advantage of much faster convergence than finite-difference or 
finite-volume methods (i.e., exhibiting exponential rather than algebraic convergence with mesh re- 
finement). These computations were performed with the same spectral computer code used to obtain 
the results presented in the paper by Zang et al. [70]. The conditions of the problem considered in 
the present work are for one of the cases included in the paper of Zang et al. In the DNS computer 
code, the spatial discretization is Fourier collocation with evenly spaced points in all coordinate 
directions. The code employs a time-step splitting technique that allows sound waves to be inte- 
grated exactly in time [71], The other terms are advanced with a third-order Runge-Kutta explicit 
time-stepping scheme. An isotropic truncation is applied at the end of each Runge-Kutta stage. This 
ensures that only the resolvable modes (i.e., k < N/2 where k is the Fourier wavenumber vector) 
are nonzero. The semi-implicit scheme provides important computational advantages. For example. 
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at small fluctuating Mach numbers it can allow the use of time steps several times larger than those 
possible with a fully explicit scheme. 


4.1.2 Initial Solution 


The method for creating the initial turbulent flowfield condition is taken directly from Erlebacher 
et al. [71,72], and is similar to a method introduced by Passot and Pouquet [73]. Random fields in 
physical space are generated for density, velocity, and temperature (p, u, and T). Subsequently the 
velocity field is decomposed in Fourier space into incompressible and compressible components: 

U 0 (x) = Uq(x) + Uq (x) (63) 

The Uq ( x ) has mean value of zero by construction, and Uq(x) is divergence free. These components 
are transformed back to physical space and re-scaled in order to impose a prescribed autocorrelation 
spectrum. The spectrum 

E(k) = (k + l) 4 exp ^— 2 - K ^ , k > 1 (64) 

is imposed on each of the variables p, ftg, Uq , and T (for k = 0 E(k) = 0). Also, the initial value 
for the ratio of the compressible kinetic energy to the total kinetic energy, which is often called \, is 
set to 0.2 (see Zang et al. [70]). Note that Eq. (64) is slightly different than that reported in Zang et 
al. [70]. The reference pressure is chosen (based on a reference density, ~p re f, and reference velocity, 
U re f ) to be p re f = PrefU'^ej, so the initial nondimensional pressure is: 


__ PT 
P 7 M? ef 

The reference Mach number 


(65) 



is taken to be M re f = 0.2 for all the current simulations. The Reynolds number, which is defined as 
Re = (p re fU re fL re f) / p. re f, is taken to be 250, and the Prandtl number and ratio of specific heats 
have the standard values of Pr = 0.72 and 7 = 1.4, respectively. The reference length L re f in this 
case is unity. 

The Taylor microscale Reynolds number is defined by: 

Rex = ^XUrms (67) 

w 

where the angle brackets indicate an average over all cells in the domain, and 

A = - (Au + Xy + Aa;) (68) 



For all the cases herein, the initial Re\ is approximately 125. The x-, y-, or ^-direction Taylor 
microscale Reynolds numbers can also be found. For example, the ^-direction Taylor microscale 
Reynolds number is defined by: 

Re\, x = ^\ x xfW) ( 71 ) 

The initial Re\ }X = Re\ }V = Re\, z « 75. 
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4.1.3 Turbulence Statistics 


This section defines some of the turbulence statistics. The total turbulent kinetic energy is: 

k = \ E (“ 2 + 52 + ™ 2 ) (72 ) 

cells 

The total turbulent dissipation is: 

ceZis 


The skewness factor of the velocity derivative, which supplies a measure of the nonlinear vortex 
stretching, is given by: 


S k = 


1 

3 


((I) 3 ) 

K(ff ) 2 >] 3/2 


((|) 3 ) 

K(ff ) 2 )] 3/2 


((f) 3 ) 
[((f ) 2 )] 3/2 


(74) 


Some discussion about the skewness factor and its relationship to the numerical scheme is given in 
Appendix B. 

The enstrophy is defined by: 

enstrophy = i(|V x u\ 2 ) (75) 


The turbulence spectra is computed as follows. Given (at an instant in time) the velocities every- 
where in the field, a 3-D fast Fourier transform (FFT) is performed on each velocity component. The 
energy E'(k) within each energy “bin” is computed by summing over the squares of all the Fourier 
indices that contribute to that energy level. Then, the result is normalized via: 


E{k) 


E\k)/ / E\ K )dK 


(76) 


4.2 Results with Upwinding and Central Differencing 

Table 1 summarizes the computer runs performed using CFL3D. Both upwinding and central differ- 
encing methods were employed, and the effects of grid size, time step, and Smagorinsky filter con- 
stant c s were explored. The time step. At* , was nondimensionalized using At* = AtU re f / L re f, 
where U re f and L re f are the reference velocity scale and reference length scale, respectively. This 
time scale is a measure of “eddy turnover time” in terms of mean flow quantities. Most of the runs 
were performed at a nondimensional time step of At* = 0.002. As can be seen in the table, one run 
was performed at a time step that was one fourth the nominal time step. Although not shown, results 
using the two different time steps were essentially identical, so the larger time step. At* = 0.002, 
was deemed to be fine enough to accurately capture the temporal development with the current 
scheme in CFL3D. 


4.2.1 Implicit LES (No Subgrid Model) 

Results using CFL3D with the upwind-biased k = 1/3 scheme are shown in Figs. 3(a)-(c) for four 
different grid cell sizes. Note that the k and e were normalized by their initial values. All cases 
except the 192 3 were carried out to total time of t* = 10. The main point to be made is that on the 
coarsest grid (32 3 ) the upwind result is very poor, particularly in the fact that it shows no increase 
at all in e for 0 < t* <~ 1.5. As will be shown later, e should increase to a level of nearly 2 (near 
t* = 1.3) before dropping off. This rise in e occurs during the energy redistribution phase of the 
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Table 1 . Summary of CFL3D Cases 


K 

Grid size (# cells) 

At* 

c s 

1/3 (upwind) 

“32 s 

0.002 

0 

1/3 (upwind) 

64 3 

0.002 

0 

1/3 (upwind) 

128 3 

0.002 

0 

1/3 (upwind) 

192 3 

0.002 

0 

1 (central) 

32 3 

0.002 

0 

1 (central) 

64 3 

0.002 

0 

1 (central) 

128 3 

0.002 

0 

1 (central) 

192 3 

0.002 

0 

1 (central) 

32 3 

0.002 

0.2 

1 (central) 

64 3 

0.002 

0.2 

1 (central) 

128 3 

0.002 

0.2 

1 (central) 

32 3 

0.002 

0.1 

1 (central) 

32 3 

0.002 

0.4 

1 (central) 

32 3 

0.0005 

0.2 


simulation. The skewness using upwind with no model ends up between about —0.2 and —0.5. For 
incompressible flows the experimental value for isotropic turbulence is around —0.4 (Lesieur [74]), 
and the DNS value at Re\ = 40 is —0.5 (Erlebacher [71]). 

Isocontours of vorticity are shown at t* = 5 on three of the grids in Fig. 4. (Contours on the 
192 3 grid are not shown; visually there are more small scale structures than for the 128 3 grid, as 
expected. The entire cube, 27 t on a side, is shown here along with the subdivisions used for parallel 
processing.) However, in the literature, the so-called Q-invariant or Q-criterion is often preferred for 
visualizing LES results. See, for example, Piomelli et al. [75] and Slomski and Chang [76]. The Q 
variable is the second invariant of the velocity gradient tensor: 


Q 


1 

2 


^ 2 j ^ 


ij 


Rij Rij ) 


1 / dut duj \ 

2 \cfEj dxi ) 


(77) 


The Q-invariant is useful because it tends to highlight the smaller-scale structures. Figure 5 shows 
the Q-invariant plotted for the upwind results. It is evident that using upwind on the 32 3 grid there 
is little eddy content, whereas on finer grids there are additional smaller scale structures present. 

Results using CFL3D with second-order central differencing (k = 1) and no SGS model (c s = 
0) are shown in Figs. 6(a) -(c) for four different grid sizes. These results can be thought of as 
under-resolved DNS computations because there is no numerical or SGS model dissipation, and the 
grids are too coarse for a second-order finite-volume scheme to resolve down to the Kolmogorov 
dissipation length scale. Comparing Fig. 6(b) with Fig. 3(b), the central differencing yielded better 
results for temporal development of e on a given grid size (the peak level was closer to the expected 
value of approximately 2). This improvement is to be expected because upwinding is inherently 
more dissipative. For example, using central differencing, the peak normalized e reached nearly 1.4 
on the 64 3 grid, which was almost as good as the upwind result on the 128 3 grid. However, note 
that the skewness tended to be smaller in magnitude using central differencing compared to upwind. 
The reason for this underprediction is not known. 

Isocontours of the Q-invariant using central differencing are shown at t* = 5 on three of the 
grids in Fig. 7. Comparing to Fig. 5, it is easy to see that central differencing yielded greater eddy 
content than upwind differencing on a given grid size. 

Figure 8 shows energy spectra for eight different results at t* = 5. The theoretical slope of 
k -5/s f rom Kolmogorov theory (known as the Kolmogorov —5/3 law) is also indicated in the plot. 
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In the inertial subrange where the energy cascade transfers energy from large to small scales, the 
spectrum should follow this slope. However, very high Reynolds numbers are required in order to 
see evidence of the inertial subrange. The current Reynolds number of 250 is far too low for a region 
with constant slope to emerge. See Appendix C for additional discussion. Nonetheless, in this and 
subsequent plots of the energy spectrum, the « -5 / 3 line will continue to be shown to help provide a 
reference for the spectrum. In all cases shown in Fig. 8, results with central difference were better 
than for upwind (less rapid drop-off in energy with increasing k). The greater decay in the upwind 
result is a consequence of the increased numerical dissipation with the upwind differencing, which 
results in greater smearing of the small scale structures associated with the higher wavenumbers. 
However, the finer the grid, the closer the central and upwind results were to each other. Figure 9 
shows the progression of the energy spectrum with eddy turnover time for central differencing (no 
SGS model) on the 128 3 grid. 

4.2.2 Explicit LES with Smagorinsky Model 

Central differencing was also used in conjunction with the Smagorinsky model, using CFL3D. 
Figs. 10(a) -(c) show the statistical information over time. Generally, the higher c s , the worse 
the result (on the 32 3 grid). The best result was achieved using no model at all (c s = 0). This 
conclusion is also reflected in the spectrum results at t* = 5, which is shown in Fig. 11. The higher 
the c s coefficient, the more the energy dropped off at the higher wavenumbers. 

On finer grids, the effect of c s was not as great, as seen in Figs. 12(a) -(c) and Fig. 13. For 
example, on the 128 3 grid, the energy spectrum result using c s = 0.2 is only slightly worse than the 
result with c s = 0. Isocontours of the Q-invariant using central differencing with the Smagorinsky 
model and c s = 0.2 are shown at 0 5 in Fig. 14. 

Upwind was not run in conjunction with the Smagorinsky model (c s > 0), because the combi- 
nation would yield even more dissipative results than upwind ILES ( c s = 0). 


4.3 Implicit LES Results with WENO and WENO-L 

The focus in this section is on the effect of higher order accuracy on the turbulence statistics. Fifth- 
order, seventh-order, and ninth-order schemes are all considered. To ensure that the best possible 
results are obtained with the higher order approximations, we have switched off the nonlinear part 
of the WENO scheme for the initial set of computations. The linear scheme is designated WENO-L. 
In the second part of this section we investigate the nonlinear scheme, which is simply designated 
as WENO. The two ways of computing the nonlinear weights, u>k, and smoothness measure, (3k, de- 
scribed in Section 3.2.1, are considered. In addition, the effect of the parameter e w is also evaluated 
and discussed. It is demonstrated by numerical testing that nearly the same solutions can be obtained 
with the nonlinear WENO scheme as for the WENO-L scheme if a sufficiently large value for e w 
in the non-oscillatory weights is chosen. For all high-order schemes the effect of mesh resolution 
on the computed solutions is examined. The grid density varied from 32 3 to 192 3 ; these are the 
same grids used for all results in the previous section. Results obtained with the WENO schemes 
are also compared with the second-order central (c s = 0 and c s = 2) and third-order upwind biased 
schemes. In addition, they are compared with those from the spectral scheme of Zang et al. [70]. For 
all simulations a At* = 0.002 was used. Reducing the time step by a factor of four, which reduces 
temporal error by a factor of 64, produced essentially the same numerical results. All of the WENO 
and WENO-L calculations were performed without an explicit SGS model. 

4.3.1 Linear WENO Scheme 

Figures 15-17 show the turbulence quantities computed using the WENO-L scheme with different 
orders of spatial accuracy. Comparing the decay in the turbulent kinetic energy of the dissipative 
higher order schemes with that of the nondissipative central scheme, which was shown in Fig. 6(a), 
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Table 2. Normalized Dissipation Peak 


Scheme 

Order 

Grid size (# cells) 

Dissipation Peak 

Central 

2 

32 s 

1.10 

Central 

2 

64 3 

1.40 

Central 

2 

128 3 

1.64 

Central 

2 

192 3 

1.75 

Upwind 

3 

32 3 

1.00 

Upwind 

3 

64 3 

1.10 

Upwind 

3 

128 3 

1.46 

Upwind 

3 

192 3 

1.65 

WENO-L 

5 

32 3 

1.00 

WENO-L 

5 

64 3 

1.37 

WENO-L 

5 

128 3 

1.75 

WENO-L 

5 

192 3 

1.87 

WENO-L 

7 

32 3 

1.08 

WENO-L 

7 

64 3 

1.56 

WENO-L 

7 

128 3 

1.84 

WENO-L 

7 

192 3 

1.91 

WENO-L 

9 

32 3 

1.17 

WENO-L 

9 

64 3 

1.69 

WENO-L 

9 

128 3 

1.96 

WENO-L 

9 

192 3 

2.02 


we observe that the dissipative schemes approach the solution on the 192 3 grid from below while 
the nondissipative scheme approaches from above, as expected. Moreover, it is evident that the 
fifth-order scheme is highly dissipative on the coarser meshes, especially on the 32 3 mesh. In (b) 
of Figs. 15-17, we see that the peak value of the normalized dissipation, which occurs during the 
energy redistribution phase of the simulation, increases both with grid refinement and with order of 
accuracy. The peak values for the different approximation orders and the range of mesh densities 
are given in Table 2, along with results from the earlier central (c s = 0) and upwind computations. 
These results are also shown in Fig. 18. The largest energy dissipation (e) peaks, obtained with the 
seventh-order and ninth-order schemes on mesh densities 128 3 and 192 3 , are roughly in the range 
of 1.8 to 2.0 and occur at t* « 1.30. On the 128 3 grid, the peak value in £ is about 20% larger 
than the one computed with second-order central differencing. Like the third-order upwind biased 
scheme, the fifth-order WENO-L scheme does not exhibit any dissipation peak on the 32 3 grid, even 
though its dissipation is somewhat smaller (as seen by comparing the turbulent kinetic energy levels 
in Figs. 3(a) and 15(a)). Even the seventh-order scheme only has a very small peak on this grid, and 
this peak is about the same size as that for the central difference scheme. The peak in dissipation for 
the fifth-order scheme exceeds 1.8 on the 192 3 grid, whereas for the seventh-order and ninth-order 
schemes it is greater than 1.8 on the 128 3 grid. These results indicate that the resolution on the 32 3 
grid is not adequate for any of the schemes. The numerical viscosity is not sufficiently small on such 
a coarse grid. 

We now consider the skewness factor (Sk) plots in Figs. 15-17. On the finest mesh the quasi- 
equilibrium value is between -0.5 and -0.6 for the different approximation orders of the WENO-L 
scheme. It is nearly -0.6 using the ninth-order scheme on the 192 3 grid. This is about the same Sk 
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obtained with the third-order upwind biased scheme on the finest mesh. 

In Fig. 19, the energy spectra are shown for computations with different orders of accuracy and 
grid sizes. The fifth-order results are comparable to those obtained with the third-order upwind 
biased scheme, although the fifth-order results are more dissipative at the highest frequencies. In 
comparing the spectra for the higher order schemes with the fifth-order scheme, we observe roughly 
an order of magnitude less drop in the energy decay at the the cut-off wavenumber (N/2, where N is 
the number of mesh points), reflecting the reduced dissipation of the small scales. Figure 20 clearly 
shows how the high-frequency content improves with increasing mesh density for the ninth-order 
scheme. In Fig. 21, we see that the spectra obtained on the 128 3 grid with the ninth-order WENO 
and the second-order central ( c s = 0) schemes are nearly the same. The fifth-order WENO and 
third-order upwind spectra are also similar to each other 

Figure 22 displays the isocontours of the nondimensional Q-invariant for the ninth-order WENO- 
L scheme computations for various grid sizes. Comparing with the corresponding contours obtained 
with second-order central differencing ( c s = 0) in Fig. 7, we note that in these plots the two schemes 
appear to have similar eddy content on the 128 3 grid. However, comparing with the explicit LES 
scheme with central differencing and c s = 0.2 (see Fig. 14), we observe that there is a significant 
increase in the eddy content for ninth-order over the range of meshes. 

4.3.2 Nonlinear WENO Scheme 

Appropriate determination of the nonlinear weights is essential since these weights allow non- 
oscillatory shock capturing and determine the rate at which the nonlinear WENO scheme returns 
to a linear form. We first evaluate the two methods for computing the nonlinear weights, u>k, and 
smoothness measures, (3k- Figure 23 compares the two methods discussed in Section 3.2.1. In the 
figure, ot\ refers to the method of Jiang and Shu (Eqs. (37) and (38)), and a 2 refers to the method 
of Borges et al., Eq. (48)). For each method, two values of e w , the parameter for ensuring that 
/ 3k is bounded, are considered. The computational results are for the fifth-order and seventh-order 
schemes on the 128 3 grid. We observe in the figure a slight improvement in the prediction of the 
peak value of e with the method of Borges et al. (o:->)- For the larger s w = 10 -2 , the energy spectrum 
shows that the WENO scheme is dissipative at the highest frequencies, whereas there is some energy 
buildup at the highest frequencies for the smaller e w = Iff 6 . This behavior is most prominent for 
the fifth-order scheme. 

In Fig. 24, the method of Borges et al. is considered, and the effect of varying e w on the tur- 
bulence statistics using a 64 3 grid is shown. The WENO and WENO-L results are nearly the same 
when e w = 10 -1 . However, such a large value could easily adversely affect the nonlinear behavior 
(i.e., shock capturing properties) of the WENO scheme. Shen and Zha [57, 58] have demonstrated 
for transonic fluid dynamics problems with shock waves that a value of e w = 10“ 2 is satisfactory. 
It still remains, however, to show that this value of s w is suitable for a wide range of flow problems 
with various levels of complexities and Mach numbers. To reduce the possibility of energy buildup 
at the highest frequencies that can result in instability, the method of Borges et al. with s w = 10 -2 
will be used for all subsequent results. 

Figure 25 compares the dissipation computed with the WENO-L and WENO schemes for dif- 
ferent grids and orders of discrete approximation. In general, the linear and nonlinear results are 
quite close, with the greatest difference occurring on the coarsest (32 3 ) grid. A similar comparison 
is shown in Fig. 26 for the energy spectrum. There is no discernible difference between the linear 
and nonlinear results. 

We now consider how the ninth-order WENO scheme compares with other previously presented 
LES schemes. Figure 27 displays the turbulence statistics for the ninth-order WENO ILES scheme, 
second-order central scheme with the Smagorinsky SGS model (c s = 0.2), and the third-order 
upwind biased ILES scheme. Results are given for the 64 3 and 128 3 grids. The largest differences 
in the turbulent kinetic energy (k) decay for the two grids occur with the upwind scheme. With the 
WENO scheme the k decay for the 64 3 grid closely tracks that for the 128 3 grid. When comparing 
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the dissipation (e) peaks, we see that the peak on the 128 3 grid with WENO scheme is significantly 
larger than the peaks for the lower-order LES schemes. Furthermore, the e peak on the 64 3 grid 
with the WENO scheme is approximately 8% and 15% higher than the peaks on the 128 3 grid with 
the central and upwind schemes, respectively. These results indicate that more than an order of 
magnitude increase in the grid density would be required to produce similar e peaks with the lower 
order LES schemes. In Fig. 27(c) a comparison of the LES schemes is given for the energy spectrum 
at t* = 5. The WENO result for the 64 3 grid follows the fine grid energy spectrum until between 
k = 15 and n = 20, where it begins to show a more rapid decrease in the energy with increasing k. 
As expected, the lower order LES schemes exhibit a much earlier departure (roughly k = 5) from 
the fine grid result with WENO. 

A comparison of the turbulence statistics at t* =2 for the ninth-order WENO ILES scheme, 
second-order central scheme (c s = 0), and a spectral scheme is presented in Fig. 28. In the absence 
of numerical dissipation (either through discretization or filtering) and an explicit SGS model, the 
central scheme relies on the dissipative effect of the physical viscous terms for stability. Thus, the 
computation with the central scheme is essentially an under-resolved DNS. One could view this 
solution as representing a bound on the best possible result with a second-order scheme for an 
explicit SGS model (i.e., in the limit as c s — > 0). The results for the spectral scheme are from DNS 
computations with the computer code of Zang et al. [70]. In Fig. 28a, we observe a similar decay 
rate of the turbulent kinetic energy for all schemes. Figure 28(b) shows that the peak dissipation 
for the ninth-order WENO scheme on the 192 3 grid is approximately 4% high compared with the 
DNS, whereas the second-order central scheme is approximately 10% low. Although not shown, a 
simulation on a 384 3 grid was performed to determine if there is any further increase in the peak 
dissipation with the WENO scheme. Essentially the same peak value was obtained on the finer grid. 
In Fig. 28c, we observe fairly good agreement between the energy spectra computed with the WENO 
and spectral schemes, except at the highest frequencies where the WENO scheme exhibits a faster 
decay in k due to dissipative effects. There is also some difference at the lowest frequencies. The 
source of this difference appeared to be a consequence of variations in the methods used for the two 
schemes to determine the spectrum from the simulation data. The second-order central scheme on 
the 192 3 grid exhibits less energy content at the highest wave numbers than the ninth-order WENO 
scheme, and thus, has poorer agreement with the DNS spectrum. 

The computational cost of the WENO scheme with increasing approximation order is espe- 
cially important due to the high operation count of the scheme. A comparison of the results for the 
ninth-order and fifth-order schemes indicates that the ninth-order scheme requires half as many grid 
points in each coordinate direction as the fifth-order scheme to obtain comparable resolution. We 
compared the computational (CPU) times of these schemes for the isotropic turbulence decay sim- 
ulation using the 64 3 grid for the ninth-order scheme and the 128 3 grid for the fifth-order scheme. 
All of the simulations were performed with parallel processing. To have approximately the same 
parallel efficiency for both schemes, we used more processors for the fifth-order scheme on the 
128 3 grid than for the ninth-order scheme on the 64 3 grid (i.e., 64 and 16 processors, respectively). 
This comparison showed that the ninth-order scheme required about 17% less computing time than 
the fifth-order scheme for comparable resolution. We should point out that Shi et al. [77] showed 
that for problems with discontinuities and complex solution structures (e.g., double Mach reflection 
and Rayleigh-Taylor instability problems) that a ninth-order WENO scheme required half as many 
mesh points in each coordinate direction to obtain comparable resolution (i.e., similar numerical 
viscosity) as the fifth-order scheme. Furthermore, they showed that with parallel processing that the 
ninth-order scheme could reduce the CPU time between 66% and 80% with their implementation 
of the WENO scheme. Their results suggest that it may be possible to significantly improve the 
present implementation of the WENO scheme, especially since no attempt was made to optimize it. 
In addition, their results emphasize the capability of the high-order WENO schemes to capture small 
scale structures of complicated flows that would require much finer grids with lower order schemes. 
This is especially important in the simulation of turbulence. 
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4.3.3 Future Work: Improved WENO Scheme 

As we have seen in the previous section, the conventional fifth-order WENO scheme of Shu [36] 
with either the weights of Jiang and Shu [53] or Borges et al. [55] is too dissipative. In addition, 
there is no energy stability proof (guaranteeing all eigenvalues lie in the left half of the complex 
plane) for the conventional WENO scheme. To improve this situation Yamaleev and Carpenter [78] 
- [80] have introduced energy stable modifications for existing WENO schemes. Such modified 
schemes are called energy stable WENO (ESWENO) schemes. By modifying the nonlinear weights 
of the conventional WENO schemes, Yamaleev and Carpenter achieved a reduction of the nonlinear 
dissipation in the neighborhood of discontinuities, while retaining an energy stability of the spa- 
tial discretization. With their systematic methodology (see [78] and [80]) they have constructed 
ESWENO schemes for both upwind biased and centered finite differences. Thus, stable higher 
order central difference schemes can be used in smooth regions of the flow field, eliminating the 
dissipation in these regions. 

An essential constraint in minimizing numerical errors produced by dissipation is to enhance 
the stability of the spatial discretization. For a given order of approximation, one approach that 
can lead to not only reduction of numerical dissipation but also minimization of aliasing error is to 
use a skew-symmetric form of the convective operator (e.g., see Kravchenko and Moin [81]). This 
form of the convective operator is defined by a linear combination of the conservative (divergence) 
and non-conservative operators using a central difference discretization. Such a discrete operator 
is conservative (i.e., can be expressed in terms of telescopic fluxes), and in the absence of external 
forces and dissipation, it produces a kinetic -energy preserving scheme. With kinetic energy con- 
servation, one can show that the scheme is energy stable, which means that without dissipation the 
energy of the discrete system cannot increase (with dissipation the energy decreases). Moreover, 
in conserving kinetic energy there is no spurious production or dissipation of kinetic energy due to 
the discretization of the nonlinear convection. Some examples of papers with schemes based on a 
skew-symmetric operator are those by Lilly [82], Ducros et al. [83], Verstappen and Veldman [84], 
Feiereisen et al. [85], and Kok [86], 

While the ESWENO formulation and alternative convective operators can significantly reduce 
numerical dissipation, and thus, improve the representation of the entire energy spectrum, the dis- 
persion errors still need to be controlled. As Ghosal [87] has demonstrated, even a centered scheme 
with eighth-order accuracy requires some technique such as prefiltering to effectively remove such 
errors so that the modeling errors dominate. One possible approach for minimizing the numerical 
dispersion errors at high wavenumbers is to construct a central difference approximation with the 
dispersion-relation preserving (DRP) scheme of Tam and Webb [88]. The DRP scheme uses a larger 
stencil than necessary for a given order of approximation to create an additional degree of freedom 
to allow minimization of dispersion. 

In future work, we plan to consider both an ESWENO scheme and the skew-symmetric form 
of the convective operator. For the viscous terms we will replace the current discretization with the 
compact approximation of Kamakoti and Pantano [59], which ensures damping of the 7r-mode. We 
will also consider minimization of dispersion error. 


5 Concluding Remarks 

Numerical simulations of decaying homogeneous isotropic turbulence have been performed with 
both low-order and high-order spatial discretization schemes. Second-order central and third-order 
upwind biased differencing have been used for the low-order schemes. The WENO framework 
has been applied to obtain fifth-order, seventh-order, and ninth-order discretizations. Large eddy 
simulation results have been presented for the various schemes on four grids, ranging from 32 3 
to 192 3 . Turbulence statistics, such as kinetic energy, dissipation, and skewness, along with the 
energy spectra from simulations of the decaying turbulence problem have been used to examine and 
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compare the schemes. 

Through the simulation results for the fifth-order WENO scheme and the second-order central 
scheme, with or without the Smagorinsky SGS model (c s < 0.2), we have shown that the fifth- 
order WENO is too dissipative on the coarsest grid considered (32 3 ). While the fifth-order WENO 
was still too dissipative on the on the 64 3 grid when compared with the central scheme without a 
SGS model, this was not the case when compared with the central scheme with a SGS model and 
c s > 0.2. On the coarser grids, the fifth-order WENO scheme was somewhat less dissipative than 
the third-order upwind biased scheme. Comparisons of the fifth-order WENO scheme with the lower 
order LES schemes have indicated that the fifth-order WENO has reasonable levels of dissipation 
on the finer grids (128 3 and 192 3 ). 

With the seventh-order and ninth-order WENO schemes, we have observed a significant im- 
provement in accuracy relative to the lower order schemes, as revealed by the computed peak in the 
energy dissipation variation with eddy turnover time. The ninth-order WENO scheme has shown the 
best overall results when comparing with the other LES schemes (i.e., second-order central scheme 
with Smagorinsky SGS model and third-order upwind biased scheme). Moreover, it was indicated 
(comparing results from 64 3 and 128 3 grids) that the other LES schemes would require more than 
an order of magnitude increase in mesh density to produce a similar dissipation peak to that ob- 
tained with the ninth-order WENO Scheme. The ninth-order WENO scheme, which is an ILES 
scheme, has also been compared with a spectral scheme for DNS. Similar results were obtained for 
the turbulent kinetic energy and turbulent dissipation. Generally, the energy spectrum of the WENO 
scheme closely tracked that of the DNS spectral scheme, except at the highest frequencies, where 
the WENO scheme exhibited some numerical dissipation effects. 

The computational cost of the WENO scheme is high. However, current simulation results have 
provided indications that for the ninth-order scheme this additional cost relative to conventional 
lower order LES schemes (such as those based on central differencing with the Smagorinsky model 
and third-order upwind biased differencing) as well as fifth-order WENO is compensated by im- 
proved resolution. In addition, as revealed by Shi et al. [77], the ninth-order WENO scheme has 
the additional advantage of capturing small scale structures of complex flows that can include shock 
waves with significantly coarser meshes than required by much lower order schemes. 

Although we have observed improved accuracy with the upwind biased higher order schemes, 
there is still a need to reduce the numerical dissipation, as revealed by the computed energy spectra. 
Recent developments for WENO schemes such as ESWENO suggest that this can be done while 
retaining appropriate shock capturing capability. In addition, these schemes could possibly be de- 
signed so that they have low dispersion errors and preserve properties such as kinetic energy. Such 
schemes would then have sufficiently low numerical errors that the dominant errors in a simulation 
would be due to modeling, allowing a better understanding of SGS modeling effects that can lead to 
improved LES results. 
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Figure 1. Uniform computational grid for WENO scheme. 
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Figure 2. The five point WENO stencil constructed from three-point stencils So, Si, and S%. 
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Figure 3. Time history of various statistics for upwind ILES for various grid sizes; (a) normalized 
kinetic energy, (b) normalized dissipation, (c) skewness. 
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128 3 



Isocontours of vorticity = 0.08, 
colored by velocity magnitude 
t*=5 



Figure 4. Isocontours of nondimensional vorticity for upwind ILES for various grid sizes, at t* = 5. 
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128 3 



Isocontours of Q-invariant = 0.03, 
colored by velocity magnitude 
t*=5 



Figure 5. Isocontours of nondimensional Q-invariant for upwind ILES for various grid sizes, at 

t* = 5. 
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Figure 6. Time history of various statistics for central (c s = 0) for various grid sizes; (a) normalized 
kinetic energy, (b) normalized dissipation, (c) skewness. 
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128 3 



Isocontours of Q-invariant = 0.03, 
colored by velocity magnitude 
t*=5 



Figure 7. Isocontours of nondimensional Q-invariant for central (c s = 0) for various grid sizes, at 

t* = 5. 
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Figure 8. Energy spectrum for upwind ILES and central (c s = 0), at t* = 5; U indicates upwind, C 
indicates central. 
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Figure 9. Energy spectrum at three different times for central (c s = 0), 128 3 grid. 
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Figure 10. Time history of various statistics for central differencing with the Smagorinsky model on 
32 3 for various c s values; (a) normalized kinetic energy, (b) normalized dissipation, (c) skewness. 
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Figure 11. Energy spectrum for the Smagorinsky model with different values of c s , 32 3 , central 
differencing, at t* = 5. 
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Figure 12. Time history of various statistics for central differencing with the Smagorinsky model 
with different values of c s on 64 3 and 128 3 ; (a) normalized kinetic energy, (b) normalized dissipa- 
tion, (c) skewness. 
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Figure 13. Energy spectrum for the Smagorinsky model with different values of c s on two finer grid 
sizes, central differencing, at t* = 5. 
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128 3 



Isocontours of Q-invariant = 0.03, 
colored by velocity magnitude 
t*=5 



Figure 14. Isocontours of nondimensional Q-invariant for central differencing (Smagorinsky model 
with c s = 0.2) for various grid sizes, at t* = 5. 
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(a) (b) 



Figure 15. Time history of various statistics on different grid sizes using 5th order WENO-L ILES 
scheme; (a) normalized kinetic energy, (b) normalized dissipation, (c) skewness. 
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Figure 16. Time history of various statistics on different grid sizes using 7th order WENO-L ILES 
scheme; (a) normalized kinetic energy, (b) normalized dissipation, (c) skewness. 
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Figure 17. Time history of various statistics on different grid sizes using 9th order WENO-L ILES 
scheme; (a) normalized kinetic energy, (b) normalized dissipation, (c) skewness. 
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normalized peak e 



Figure 18. Normalized peak e as a function of grid spacing squared. 
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(a) 


(b) 




(c) (d) 


Figure 19. Energy spectrum for 5th-order, 7th-order, and 9th-order WENO-L ILES on four different 
grids at t* = 5; (a) 32 3 , (b) 64 3 , (c) 128 3 , (d) 192 3 . 
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Figure 20. Energy spectrum for 9th-order WENO-L ILES at t* - 5 
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Figure 21. Energy spectrum for C2 = central 2nd order, c s = 0 (CFL3D), U3 = upwind 3rd order 
ILES (CFL3D), W9 = WENO-L 9th order, and W5 = WENO-L 5th order ILES on same 128 3 grid, 
at t* = 5. 
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128 3 



Isocontours of Q-invariant = 0.03, 
colored by velocity magnitude 
t*=5 



Figure 22. Isocontours of nondimensional Q-invariant for 9th-order WENO-L for various grid sizes, 
at t* = 5. 
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(a) 



(c) 
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Figure 23. Effect of two different ways of computing nonlinear WENO weights on the normalized 
dissipation and energy spectrum (128 3 ); (a) dissipation, 5th order, (b) spectrum at t* = 2, 5th order 
(c) dissipation, 7th order (d) spectrum at t* = 2, 7th order. 
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Figure 24. Time history of statistics and energy spectrum comparison for 5th-order WENO (64 3 ) 
with s w = 0.1, e w = 0.01, e w = 0.001; (a) normalized kinetic energy, (b) normalized dissipation, 
(c) energy spectrum at t* = 1.5. 
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Figure 25. Comparison of normalized dissipation computed with WENO-L and WENO; (a) 5th 
order, b) 7th order, (c) 9th order. 
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Figure 26. Comparison of energy spectrum computed with WENO-L and WENO, t* = 2; (a) 5th 
order, b) 7th order, (c) 9th order. 
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(c) 

Figure 27. Time history of statistics and energy spectrum comparison for 9th-order WENO, 2nd- 
order central with Smagorinsky SGS model (c s = 0.2), and 3rd-order upwind biased; (a) normalized 
kinetic energy, (b) normalized dissipation, (c) energy spectrum at t* = 5. 
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(a) (b) 
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Figure 28. Time history of statistics and energy spectrum comparison for 9th-order WENO, 2nd- 
order central ( c s = 0), and spectral DNS; (a) normalized kinetic energy, (b) normalized dissipation, 
(c) energy spectrum at t* = 2. 
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Appendix A 

For the seventh-order WENO scheme, the fourth-order fluxes at a++ 1/2 are as follows: 

fi+ 1/2 = ^(-3/i-3 + 13/i-2- 23/^1 + 25/0, 

fl+ 1/2 = — 5/i_i + 13/j + 3/ i+ i), 

/f+ 1/2 = ^( — /»-i + 7/, + 7/j+i — 0+ 2 ), 

/f+1/2 = + 13/i+l — 5/,: + 2 + 0+ 3 ), 

The corresponding ideal weights are as follows: 

w 0 = 1/35, wi = 12/35, w 2 = 18/35, w 3 = 4/35. (Al) 

Then, substituting the // +1 / 2 and the ideal weights into Eq. (29), we obtain the seventh-order flux at 
x i+ 1 / 2 , which is given by 

O+ 1/2 = ^(-3/i-s + 25/i_ 2 - 101/i-! + 319/j + 214/ i+1 - 38 f i+2 + 4/ i+3 ). (A2) 

Using Eq. (38) and the coefficients of the third-degree interpolating polynomials, the smoothness 
indicators for the nonlinear weights of the seventh-order scheme are found to be 

A = / i _ 3 (547/ i _ 3 -3882/ < _ 2 + 4642/^!- 1854/0 

- /i- 2 (7042/j_ 2 - 17246/j_! + 7042/,) 

+ /, ,(11003./) q - 9042/0 + 2107// 

A = /i- 2 (267/j_ 2 — 1642/j_! + 16020 — 494/ i+1 ) 

- /i- 1 (2843/j_ 1 - 5966.0 + 1922/ i+1 ) 

+ 0(3443/, - 2522 f i+1 ) + 547 /? +1 

(A3) 

0 2 = 0-1(5470.! - 25220 + 1922/, +1 - 494/ i+2 ) 

- 0(34430 - 5966/ i+1 + 1602/ i+2 ) 

+ 0 +1 (28430+! - 16420+a) + 267/^ 2 

A = 0(21070 - 9402/ i+1 + 7042/ i+2 - 1854/ i+3 ) 

- 0+i (110030+1 - 172460+ 2 + 4642/ i+3 ) 

+ O+ 2 (7043O+ 2 - 38820+ 3 ) + 547/^ 3 

and Ofc = Ofe/240 for k = 0, . . . , 3. 

For the ninth-order WENO scheme, the fifth-order fluxes at a+ +1 / 2 are as follows: 
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fi+1/2 = ^(12/i-4 - 63/i_3 + 137/i_2 - 163/i_! + 137/*), 

£+1/2 = ^j(— 3/»-3 + 17/»-2 — 43/j_! + 77/, + 12/ i+1 ), 

/f+1/2 = ^(2/i-2 - 13/i-i + 47/i + 27/ i+1 - 3/ i+2 ), (A4) 

/f+i/2 = ^( — 3/i_! + 27/j + 47/j +1 — 13/»+2 + 2/j +3 ), 

Zti/2 = ^(12/i + 77/j+i. - 43/ i+2 + 17/j+ 3 - 3/ i+4 ), 

The corresponding ideal weights are as follows: 

w 0 = 1/126, wi = 10/63, w 2 = 10/21, w 3 = 20/63, w 4 = 5/126. (A5) 

Then, substituting the f^ +1 / 2 and the ideal weights into Eq. (29), we obtain the ninth-order flux at 
Xi + 1 / 2 , which is given by 

fi+i/2 = ^ [12/<- 4 - 123/, _ 3 + 597/j_ 2 - 1923/^, + 5637/, (A6) 

' 7560 

+ 4125/ i+1 - 915/, +2 + 165/j +3 - f i+i ] 

Using Eq. (38) and the coefficients of the fourth-degree interpolating polynomials, the smoothness 
indicators for the nonlinear weights of the ninth-order scheme are found to be 

A = /i_ 4 (22658/i_ 4 - 20850/, _3 + 364863/, _ 2 - 288007/, _i + 86329/,) 
/,_ 3 (482963/,_3 - 1704396 /,_ 2 + 1358458/, _ 4 - 411487/,) 

+ /,_ 2 (1521393/,_ 2 - 2462076/, _i + 758823/,) 

+ /,_i(1020563/,_i - 649501/,) + 107918/, 2 

A = /,- 3 (6908/,_3 - 60871/ i _ 2 + 99213/,_ 1 - 70237/, + 18079/,+i) 

- /,_ 2 (138563/,_ 2 - 464976/,_i + 337018 /, - 88297 /, +1 ) 

+ /,_!(406293/,_! - 611976/i + 165153/, +1 ) 

+ /, (242723/, - 140251/ i+ i) + 22658/ 2 +1 

A = /,_ 2 (6908/,_ 2 - 51001/i-i + 67923/, - 38947/, +i + 8209 f i+2 ) (A7) 

- /,-i(104963/,_! - 299076/; + 179098/, +1 - 38947/, +2 ) 

+ /,(231153/, - 299076/, +i + 67923/, +2 ) 

+ /,+i(104963/,+i - 51001/i+a) + 6908/^ 2 

= /i_i(22658/i_i - 140251/, + 165153/,+i - 88297/ i+2 + 18079/, +3 ) 

- /i (242723/,; - 611976/, +1 + 337018/, +2 - 70237 / i+3 ) 

+ /,+i (406293/, +i - 464976/ i+2 + 99213/, +3 ) 

+ /, +2 (138563/, +2 - 60871/, +3 ) + 6908 f? +3 

= /,(107918/, - 649501/, +1 + 758823/, +2 - 411487/ i+3 + 86329 /, +4 ) 

- /,+i (1020563/, +i - 2462076/, +2 + 1358458/, +3 - 233007/ i+4 ) 

+ /, +2 (1521393/, +2 - 1704396/, +3 + 364863 /, +4 ) 

+ /, +3 (482963/, +3 - 208501/, +4 ) + 22658/, 2 +4 

and /3k = /3 fc /5040 for k = 0, . . . , 4. 
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Table Al. The constants c r j for 8th order. 


r 

3=0 

3 = 1 

3 = 2 

3 = 3 

3= 4 

3 = 5 

3= 6 

3 = 7 

-i 

761/280 

-1479/280 

2441/280 

-8357/840 

6343/840 

-613/168 

171/168 

-1/8 

0 

1/8 

481/280 

-499/280 

481/280 

-1007/840 

463/840 

-25/168 

1/56 

1 

-1/56 

15/56 

341/280 

-219/280 

131/280 

-167/840 

43/840 

-1/168 

2 

1/168 

-11/168 

73/168 

743/840 

-307/840 

113/840 

-9/280 

1/280 

3 

-1/280 

29/840 

-139/840 

533/840 

533/840 

-139/840 

29/840 

-1/280 

4 

1/280 

-9/280 

113/840 

-307/840 

743/840 

73/168 

-11/168 

1/168 

5 

-1/168 

43/840 

-167/840 

131/280 

-219/280 

341/280 

15/56 

-1/56 

6 

1/56 

-25/168 

463/840 

-1007/840 

481/280 

-499/280 

481/280 

1/8 

7 

-1/8 

171/168 

-613/168 

6343/840 

-8357/840 

2441/280 

-1479/280 

761/280 


This eighth-order table is an amendment to Table 2.1 of Shu [36]. As in Table 2.1, the c rj are 
weights of interpolating polynomials; k and r denote the order of approximation and the number of 
grid points to the left of point Xi, respectively. 
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Appendix B 


Measurements generally give values for the velocity derivative skewness of Sk ~ —0.5 (note that 
sometimes the skewness is defined with opposite sign, giving Sk ~ +0.5). Whether or not the 
skewness depends on Reynolds number is still debated. Velocity derivative skewness is related to 
vortex stretching , the basic mechanism of production of small scales in turbulence is as follows: The 
spectral evolution equation for homogeneous isotropic turbulence [89] is 


dE(n, t)/dt = P(k, t) — T(k, t) — D(k, t) 


(Bl) 


where E(n,t) is the energy spectrum, P(n,t) is a production spectrum, T(n,t) is the nonlinear 
transfer, and D(k, t ) = 2 vk 2 E(k, t) is the dissipation spectrum. The second moment of this equa- 
tion is the enstrophy balance [89] 

poo pOO 

/ dn n 2 dE{n,t) /dt = / dn k 2 [P(k, t) — T(k, t) — D(k, t)] 

Jo Jo 

If, as is usual, production is concentrated at large scales, then 

pOO 

/ dKK 2 P(«,f)«0 

Jo 

Writing the energy transfer as T(n , t) = dP(n, t)/dn , where T is the energy flux, 

poo pOO 

— I dti k 2 T(k) = / dtv 
Jo Jo 

pOO 

so that — dn k 2 T(k) > 0 if the energy flux is from large to small scales. Thus, in the approxi- 

J o 

mate enstrophy balance 


pOO pOO 

/ dn k 2 8E(k , t) /dt = — / dn [k 2 T(k, t) + k 2 D(k, f)] 
Jo Jo 


(B2) 


The first term on the right side is a source term that generates enstrophy. Elementary considerations 
relate this term to vortex stretching [89]. Bl The Kolmogorov theory asserts that the small scales 
in any turbulent flow are approximately isotropic and homogeneous, thus, this special analysis for 
homogeneous isotropic turbulence is believed to apply in much greater generality. 

One of the first results of Taylor’s statistical theory of turbulence [91] is 


S k = -V30 


dn k 2 T(k) 


dn k 2 E(k) 


1 3/2 


(B3) 


Then, in view of Eq. (B2), velocity derivative skewness Sk is a dimensionless measure of the produc- 
tion of enstrophy, hence of the generation of small scales. The skewness also measures the strength 
of the nonlinearity: For example, it vanishes in a Gaussian random field. The Batchelor skewness 
relation 


pOO pOO 

/ dn n 2 T{n^t) = / dtv2uK 4: E(i^,t) 

Jo Jo 


B1 That vortex stretching is in fact positive was demonstrated experimentally by Taylor, who doubted von Karman’s sug- 
gestion that vortex stretching vanishes [90]. 


65 



permits an evaluation of the skewness in terms of integrals over the energy spectrum alone: since 



e 

v ’ 


S k = - — V30 


dn 2 isk 4 E(k) 

W 2 

dn k 2 E(k) 


3 

14 




3/2 


dn 2 vk 4 E{k) 


(B4) 


Explicit evaluation requires knowing the exact energy spectrum, but the trends can be assessed 
by assuming a model spectrum of the form 


E(k) = C K 


e 2/3 K -5/3 

eV3 K d 5/3 f( K / K d) 


K < Kd 
K> K d 


(B5) 


where K d = a(e/i' 3 ) 1 / 4 is proportional to the Kolmogorov microscale and Ck ~ 1-5 is the Kol- 
mogorov constant. Experimental and numerical data support the dissipation range form 


/(«/Kd) 



b 

exp[a(l 


K/n d )\ 


Direct interaction approximation (DIA) predicts this form with b = 2 and a « 4. [92], Assuming 
the general expression Eq. (B5), 

S k = - ^V30a w ' 3 C K Q+2 e a J dx x 4+b e~ ax ^ (B6) 

Ignoring the dissipation range contribution would result in S k ~ —ct 10 / 3 . 

The most important conclusion is that the result is dominated by the small scales because of the 
dependence on the parameter a, which determines the scale at which dissipation becomes signif- 
icant, and the spectrum in the dissipation range.® 2 However, the value of n d itself does not enter 
Eq. (B6). In this analysis, then, the skewness is independent of Reynolds number; while it remains 
uncertain whether this conclusion is supported by experimental data, we note that Reynolds number 
dependence would require a departure from the Kolmogorov spectrum. 

The value of the velocity derivative skewness is often taken as a figure of merit to validate nu- 
merical schemes; the connection of the skewness to vortex stretching and generation of small scales 
certainly makes this reasonable. The previous analysis is easily applied to the case of the “numerical 
turbulence” generated by a finite-difference scheme on a lattice with uniform isotropic spacing h. 
Following Vreman et al. [64], such numerical turbulence is characterized by the multipliers lj that 
express the action of the finite-difference operator Dj on “plane wave” exponential functions, 

= ily x 


Vreman et al. gives several examples. We cite the simplest one; the “vertex centered” operator for 
which 

sin (Kjh) 
i= h 

The previous analysis can be repeated essentially verbatim , replacing wavevector k and its cor- 
responding wavenumber magnitude k = \k\ by 1 and l = |1|. The result is a type of regularization 
of the Navier-Stokes equation in which the lattice spacing h enters as a new parameter. The spec- 
tral evolution equation remains valid, but the constant flux Kolmogorov spectrum will be modified 

DO • ..... 2 , 

Compare Smith and Reynolds [93], although their conclusion that the dissipation range should depend on e~ K instead 
of e~ K is not supported by data or theory. 
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by the regularization for scales k near the lattice spacing h: how near depends on the difference 
scheme. The value of the skewness will depend on the behavior of the numerical spectrum, and 
a value close to Sk ~ —0.5 must not be expected unless the Kolmogorov spectrum can extend to 
scales significantly larger than h. On the other hand, the sensitive dependence of the skewness to 
small scale details justifies caution in making positive or negative claims for difference schemes on 
this basis alone. 
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Appendix C 


Inertial Subrange 

A turbulent flow field contains (vortical) structures, which are generally called eddies. There are 
multiple length scales associated with the eddies of various sizes, starting with the largest eddies 
having a length scale proportional to the characteristic length of the field and terminating with the 
small eddies with a length scale proportional to the Kolmogorov length scale. Nonlinear interactions 
among eddies with a comparably large length scale (i.e., small wavenumber k) are approximately 
inviscid in nature and can result in eddies with smaller length scales (i.e., higher wavenumbers). 
One may view this process as large scale structures breaking up into smaller scale structures. 

There is a fluctuating velocity associated with each eddy; and thus, there is a certain kinetic 
energy corresponding to each eddy. Each fluctuating velocity has a frequency. The largest eddies 
correspond to low frequencies, and the smallest correspond to high frequencies. Consider a turbulent 
flow field that is statistically homogeneous with respect to time. Then, the turbulent kinetic energy 
of an eddy can be represented by a full range of wavenumbers. An ideal energy spectrum for 
homogeneous isotropic turbulence is displayed in Fig. Cl. The spectrum shows that during the 
initial increase of the wavenumber k the energy increases, due to the largest eddies in the domain, 
to a peak level, where the eddies have an integral length scale that is an 0(1) fraction of the large 
scales. Interactions between the eddies with an integral scale start an energy cascade process. The 
energy decreases as the frequencies of the spectrum increases. It decreases with a slope of —5/3, and 
this part of the energy spectrum which is characterized by the Taylor microscale is often called the 
Kolmogorov inertial subrange. To have a well defined inertial subrange there must be a sufficiently 
large Taylor microscale Reynolds number. This regime of the energy spectrum is followed by a 
region where the smallest eddies are dissipated. 

It is difficult to establish an inertial subrange with « -5 / 3 behavior in decaying isotropic turbu- 
lence. As indicated by Ishihara [45], a Taylor microscale Reynolds number Re\ > 200 is required 
to have an inertial subrange. In order to produce a significant inertial subrange and to investigate the 
physics of the small scales of isotropic turbulence, an alternative simulation approach is frequently 
considered. With this approach, a forcing function is introduced to emulate the energy transfer to 
the inertial subrange from the larger energy containing eddies. Thus, in a simulation the large scales 
are modeled and the small scales are resolved, and higher Re\ can be realized, which can result in 
a sufficiently wide inertial subrange. Due to the steady forcing the turbulence becomes statistically 
steady. 

Since the Taylor microscale Reynolds number for the case being considered is relatively low, 
and there is no forcing function, we should observe at most a very small inertial subrange. Never- 
theless, as a reference, we include the inertial subrange curve which has a -5/3 slope. This facilitates 
comparing with results from other investigators who include this curve. 
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Figure Cl. Sketch of Kolmogorov energy spectrum indicating length scales corresponding to differ- 
ent parts of the spectrum. 


69 


REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704-0188 


The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 
Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 
shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 


5c. PROGRAM ELEMENT NUMBER 


5f. WORK UNIT NUMBER 

561581.02.08.07.46.04 


11. SPONSOR/MONITOR'S REPORT 
NUMBER(S) 

NASA/TM-20 12-2 17593 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified - Unlimited 
Subject Category 02 

Availability: NASA CAS1 (443) 757-5802 

13. SUPPLEMENTARY NOTES 


14. ABSTRACT 

Numerical simulations of decaying homogeneous isotropic turbulence are performed with both low-order and high-order spatial discretization schemes. 

The turbulent Mach and Reynolds numbers for the simulations are 0.2 and 250, respectively. For the low-order schemes we use either second-order central or 
third-order upwind biased differencing. For higher order approximations we apply weighted essentially non-oscillatory (WENO) schemes, both with linear and 
nonlinear weights. There are two objectives in this preliminary effort to investigate possible schemes for large eddy simulation (LES). One 
is to explore the capability of a widely used low order computational fluid dynamics (CFD) code to perform LES computations. The other is to determine 
the effect of higher order accuracy (fifth, seventh, and ninth order) achieved with high-order upwind biased WENO-based schemes. Turbulence statistics, such 
as kinetic energy, dissipation, and skewness, along with the energy spectra from simulations of the decaying turbulence problem are used to assess and 
compare the various numerical schemes. In addition, results from the best performing schemes are compared with those from a spectral scheme. The effects of 
grid density, ranging from 32 cubed to 192 cubed, on the computations are also examined. The fifth-order WENO-based scheme is found to be too dissipative, 
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relative to the lower order LES schemes, as revealed by the computed peak in the energy dissipation and by the energy spectrum. 
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